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1  Introduction 


Breast  cancer  is  the  most  frequently  diagnosed  malignancy  among  women  in  the  United 
States  [1],  In  1996  the  American  Cancer  Society  estimates  that  184,300  women  will  be 
newly  diagnosed  with  breast  cancer  and  that  44,300  will  die  from  the  disease  [1].  Breast 
cancer  accounts  for  31%  of  all  cancers  detected  and  17%  of  all  cancer  deaths,  and  ranks  as 
the  second  leading  cause  of  death  from  cancer  among  women  in  the  United  States  [1].  Five 
year  survival  rates  are  generally  very  high  (93%)  for  breast  cancer  staged  as  being 
localized,  falling  to  72%  for  regional  disease  and  only  18%  for  distant  disease  [2].  The  early 
detection  of  breast  cancer  is  clearly  a  key  ingredient  of  any  strategy  designed  to  reduce 
breast  cancer  mortality. 

Mammography’s  role  is  the  early  detection  of  breast  cancer.  Although  more  accurate  than 
any  other  modality,  existing  techniques  for  mammography  only  find  80  to  90%  of  the  breast 
cancers.  Moreover,  in  7  to  10%  of  cases,  the  cancer  will  not  be  visible  on  the  mammogram. 
It  has  been  suggested  that  mammograms  as  normally  viewed,  display  only  about  3%  of  the 
total  information  detected.  Perception  is  a  problem  particularly  for  patients  with  dense 
fibroglandular  patterns.  The  importance  of  diagnosis  of  breast  cancer  at  an  early  stage  is 
critical  to  patient  survival.  The  general  inability  to  detect  small  tumors  and  other  salient 
features  within  mammograms  motivates  our  investigation. 

The  goal  of  this  project  is  to  develop  a  diagnostic  tool  for  radiologists  that  will  refine  the 
perception  of  mammographic  features  (including  lesions,  masses  and  calcifications)  and 
improve  the  accuracy  of  diagnosis.  Our  research  efforts  are  geared  towards  improving  the 
local  mammographic  viewing  environment  by  selectively  processing  mammograms  along 
different  orientations,  and  towards  providing  a  better  global  mammographic  viewing 
environment  by  fusing  together  locally  processed  sections  of  images.  By  improving  the 
visualization  of  breast  pathology  we  can  increase  the  chances  of  early  detection  of  breast 
cancers  (improve  quality)  while  requiring  less  time  to  evaluate  mammograms  for  most 
patients  (lower  costs). 

A  major  reason  for  poor  visualization  of  small  malignant  masses  is  the  subtle  difference  in 
x-ray  attenuation  between  normal  glandular  tissues  and  malignant  disease  [3].  This  fact 
makes  the  detection  of  small  malignancies  problematical,  especially  in  younger  women  who 
have  denser  breast  tissue.  Although  calcifications  have  high  inherent  attenuation 
properties,  their  small  size  also  results  in  a  low  subject  contrast  [4].  As  a  result,  the 
visibility  of  small  tumors,  and  any  associated  microcalcifications,  will  always  be  a  problem 
in  mammography  as  it  is  currently  performed  using  analog  film. 

We  are  investigating  a  methodology  for  accomplishing  mammographic  feature  analysis 
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through  multiscale  representations.  Wide  variety  of  feature  sizes  and  shapes  in 
mammograms  makes  single-scale  processing  methods  ineffective.  In  [5],  it  was  shown  that, 
in  the  context  of  mammography,  multiscale  image  processing  algorithms  can  outperform 
traditional  contrast  enhancement  methods  such  as  histogram  equalization  and  unsharp 
masking.  As  reported  in  [6],  an  improvement  in  feature  visualization  was  noted  for 
mammograms  processed  using  multiscale  wavelet  processing  techniques.  Furthermore,  in 
[7],  it  was  demonstrated  that  unsharp  masking  with  a  Gaussian  lowpass  filter  can  be 
formulated  as  a  special  case  of  contrast  enhancement  via  a  discrete  dyadic  wavelet 
transform.  Here,  we  pay  special  attention  to  two-dimensional  extensions  of  a  discrete 
dyadic  wavelet  transform  that  enable  efficient  directional  processing  of  mammographic 
images,  and  then  examine  the  usefulness  of  the  derived  transform  for  image  fusion  in 
mammography. 

In  the  sections  that  follow,  we  briefly  overview  the  contents  of  the  report,  list  publications, 
and  explain  notation  used  in  the  report. 

1.1  Overview  of  Contents 

Discrete  dyadic  wavelet  transform  has  been  successfully  applied  to  processing  of 
mammographic  images  [8,  5,  7,  9,  10].  Since  the  transform  uses  wavelets  that  are 
derivatives  of  central  B-spline  functions,  we  begin  Section  2  with  a  review  of  central 
B-splines’  properties  and  algorithms. 

Extensions  of  the  discrete  dyadic  wavelet  transform  are  presented  in  Section  2.2.  Since  the 
transform  is  of  paramount  importance  in  our  methodology  for  improving  the  imaging 
performance  of  mammography,  this  section  represents  the  largest  portion  of  the  report. 
First,  some  of  the  shortcomings  of  orthogonal  and  biorthogonal  wavelet  analyses  are 
presented  in  Section  2.2.1.  Lack  of  translation  and  rotation  invariance  of  these  transforms 
motivates  the  use  of  redundant  representations.  Next,  Section  2.2.2,  examines  a  dyadic 
wavelet  transform  in  one  dimension.  We  present  extensions  of  the  discrete  dyadic  wavelet 
transform  to  higher  order  derivatives  and  to  even  order  spline  functions.  To  compute  a 
discrete  transform  from  a  continuous  one,  the  discrete  computation  must  be  properly 
initialized.  We  devise  a  new  initialization  procedure  that  is  shown  to  be  more  accurate  than 
the  one  suggested  by  Mallat  and  Zhong  [11],  After  the  derivation  of  the  transform,  we  point 
out  relevant  connections  to  scale-space  filtering  and  reconstruction  from  edges  in  Section 
2.2.3.  Section  2.2.4  then  presents  a  direct  extension  of  the  one-dimensional  transform  to 
two  dimensions.  In  case  of  the  first  derivative  of  a  Gaussian,  a  rotation-invariant  transform 
can  be  obtained  by  a  tensor  product  extension  of  the  one-dimensional  transform  to  two 
dimensions.  Anisotropies  introduced  by  the  two-dimensional  second  derivative  discrete 
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dyadic  wavelet  transform  may  be  avoided  by  combining  derivatives  into  an  isotropic 
operator,  such  as  Laplacian  of  Gaussian  [7,  10].  Such  an  operator,  however,  cannot  perform 
orientation  analysis  (due  to  its  isotropic  nature).  To  lay  grounds  for  effective  directional 
operators,  the  concept  of  steerability  is  reviewed  in  Section  2.2.5.  Steerable  dyadic  wavelet 
transform  is  constructed  in  Section  2.2.6,  and  its  x-y  separable  analog  in  Section  2.2.7. 
These  transforms,  in  addition  to  being  translation-invariant,  enable  rotation-invariant 
processing  for  derivatives  of  orders  higher  than  two  as  well.  Given  that  the  developed 
transforms  are  highly  redundant,  an  efficient  implementation  is  extremely  advantageous. 
This  issue  is  addressed  in  Sections  2.2.8  and  2.2.9.  When  filtering  of  a  signal  is  performed 
by  circular  convolution,  boundary  effects  may  result.  A  standard  way  of  dealing  with  this 
problem  is  by  mirror-extending  the  input  signal  to  the  filter.  We  build  on  ideas  presented 
in  [12],  and  combine  such  an  extension  with  symmetry/antisymmetry  of  the  filters,  to 
achieve  savings  in  both  computation  time  and  memory. 

Section  2.3  deals  with  image  fusion.  First,  two  popular  transforms  for  image  fusion 
applications,  the  gradient  pyramid  and  orthogonal/biorthogonal  wavelet  transform,  are 
briefly  presented  in  Sections  2.3.1  and  2.3.2.  Comparison  of  the  two  transforms  with 
steerable  dyadic  wavelet  transform  using  subjective  and  objective  criteria  is  presented  in 
Section  2.3.3.  Steerable  dyadic  wavelet  transform  did  not  introduce  artifacts  common  for 
the  orthogonal  and  biorthogonal  wavelet  transforms,  while  outperforming  the  gradient 
pyramid  in  terms  of  sharpness  and  mathematically  defined  error  criteria. 

1.2  Publications 

Below,  we  provide  the  list  of  publications  accomplished  during  the  first  year  of  the  project. 

[1]  I.  Koren,  A.  Laine,  and  F.  Taylor,  “An  overcomplete  enhancement  of  digital 
mammograms,”  Era  of  Hope,  A  Multidisciplinary  Reporting  of  DoD  Progress ,  Washington, 
D.C.,  Oct.-Nov.  1997. 

[2]  I.  Koren  and  A.  Laine,  “A  discrete  dyadic  wavelet  transform  for  multidimensional 
feature  analysis,”  in  Time- Frequency  and  Wavelets  in  Biomedical  Signal  Engineering,  M. 
Akay,  Ed.,  IEEE  Press,  New  York,  NY,  1997,  pp.  425-449. 

[3]  S.  Schuler,  I.  Koren,  M.  Shim,  A.  Laine,  B.  Steinbach,  and  W.  Huda,  “An  interactive 
tutorial  for  contrast  enhancement  of  digital  mammograms  via  multiscale  representations  on 
the  Web,”  RSNA  EJ,  1997. 

[4]  A.  F.  Laine,  I.  Koren,  S.  Schuler,  W.  Huda,  and  B.  G.  Steinbach,  “Contrast 
enhancement  of  mammographic  features  via  multiscale  analysis,”  RSNA  82nd  Scientific 
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Assembly  and  Annual  Meeting ,  Chicago,  IL,  1996. 

1.3  Notation 

We  use  symbols  N,  Z,  and  R  for  the  sets  of  naturals,  integers,  and  reals,  respectively. 
L2(R)  and  L2(R2)  denote  the  Hilbert  spaces  of  measurable,  square-integrable  functions 
f(x)  and  f(x,y),  respectively. 

The  inner  product  of  two  functions  f(x)  G  L2{R)  and  g(x)  G  L2(R )  is  given  by 

/OO 

f(x)  g{ x)  dx. 

-OO 

The  norm  of  a  function  f(x)  G  L2(R)  is  defined  as 

ii/ii = 

The  convolution  of  functions  /( x)  G  L2{R)  and  g(x)  G  L2(R)  is  computed  as 

/OO 

f(t)  g(x-t)  dt, 

-oo 

and  the  convolution  of  two  functions  f(x,y )  G  L2(R2)  and  g(x,y)  G  L2(R2)  equals 

/OO  r  OO 

/  f(tx,ty)  g(x -tx,y -ty)dtxdty. 

-OO  J  —  OO 

The  Fourier  transform  of  a  function  f(x)  G  L2(R)  is  defined  as 

~  r°° 

/M=  /  f(x)e~^xdx, 

J — OO 

and  the  Fourier  transform  of  a  function  f(x,y )  G  L2(R2)  is  equal  to 

/>»,  w„)=  f°°  [°°  f(x,y)e-^MyyUxdy. 

i2(Z)  and  l2(Z2)  stand  for  the  spaces  of  square-summable  discrete  signals  f(n)  and 
f(nx,ny),  respectively. 

The  ^-transform  of  a  discrete  signal  f(n)  G  l2(Z)  is  defined  as 

OO 

F(z)=  E  /(«)*■"• 

Tl—  —  OO 

The  convolution  of  discrete  signals  /(n)  G  l2(Z )  and  g(n)  G  l2(Z)  is  equal  to 

00 

f*g(n)=  Y,  f(m)g(n-m), 

m~— oo 
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and  the  convolution  of  discrete  signals  f(nx,ny )  G  l2(Z2)  and  g(nx,ny)  G  l2(Z2)  is  given  by 

00  OO 

f*g(nx,ny)=  X)  / (mx,  my)  g(nx  —  mx,  ny  —  my). 

mx=— oo  my=— oo 

The  Fourier  transform  of  a  discrete  signal  /(n)  G  l2(Z)  is  equal  to  the  ^-transform 
evaluated  on  the  unit  circle 

00 

F(w)  =  E  /(»)«■*”. 

?1  =  — OO 

and  the  Fourier  transform  of  a  discrete  signal  f(nx,ny )  G  l2(Z2)  is  defined  as 

00  OO 

F(uI,»,)=  E  E  /(«.,»,) 

71  x —  —  00  Tty  —  —  00 

For  later  use,  we  define  the  following  functions: 


1.  the  unit  impulse  function 


2.  the  unit  step  function 


for  x  =  0 
otherwise, 


for  x  >  0 
for  x  <  0, 


3.  the  rectangular  function 


rect(:r) 


1  for  |a:|  <  \ 
0  for  ja;j  >  ~, 


4.  the  sine  function 


5.  the  unit  impulse  sequence 


sinc(a;) 


sin(7ra;) 
irx  ’ 


and 


for  n  =  0 
otherwise, 


where  x  G  R  and  n  G  Z. 
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2  Body 

2.1  Signal  Processing  Using  Central  B-Splines 

In  this  section,  we  briefly  review  fundamentals  of  spline  processing  needed  for  derivations  in 
Section  2.2.  First,  Section  2.1.1  presents  basic  properties  of  central  B-splines.  Next,  Section 

2.1.2  introduces  a  digital  filtering  scheme  for  B-spline  signal  processing.  Finally,  Section 

2.1.3  describes  signal  approximations  using  projections  onto  the  spline  function  spaces. 

2.1.1  Central  B-Splines:  Definition  and  Properties 

Given  real  numbers  —  oo  <  xo  <  xi  <  x2  < . . .  <  xm  <  xm+i  <  oo,  a  function  on  the  interval 
[xoj  xm+i]  is  called  a  spline  function  of  order  p  with  the  knot  (i.e.,  grid  point)  sequence 
Xi,  x2, . . .  xm,  if  it  is  (1)  a  polynomial  of  degree  p  or  less  in  each  interval  [a;,,  xi+1], 
i  =  0, 1, . .  .m,  and  (2)  continuous  in  its  derivatives  up  to  the  order  p—  1  on  the  interval 
[x0,xm+1]  (i.e.,  Cp~l[x 0,rrm+i]). 

Here,  we  will  concentrate  primarily  on  basis  splines  (B-splines),  or  more  precisely,  central 
B-splines  having  knots  at  i  G  Z  for  p  odd  and  at  i  +  |  for  p  even  [13].  Central  B-splines  of 
order  p  (with  p+1  knots)  are  defined  as 


Figure  1  shows  (3p{x )  and  their  Fourier  transforms  f3p(u)  for  p  G  {0, 1, 2, 3, 4}. 

A  family  of  functions  {/3p(x  -  m )}mez  forms  a  basis  of  Sp,  a  space  of  order  p  spline 
functions  with  knots  at  i  for  p  odd  and  at  i  +  |  for  p  even  ( i  G  Z)  [13,  14].  Except  for 
p  =  0,  the  basis  functions  {/ 3p(x  —  m)}  are  not  orthogonal. 

Let  us  list  some  properties  of  functions  (3p(x )  [15,  13]: 

1.  (3p(x)  are  nonnegative  functions  with  a  support  of  length  p+1, 

p+1  times 
/ - * - ■. 

2.  (3p(x)  =  Pq  *  p0  *  •  •  •  *  (3q(x),  where  denotes  the  convolution  operator,  or, 
equivalently,  in  the  Fourier  domain: 


3.  (3p(x)  -  J  ((^  +*)  PP- i(x  +  |)  +  -  x)  (3P-i(x  - 

4  dMzl  =  pp_^x  +  i)  _  pp_x(x  _  I). 
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Another  interesting  property  of  B-splines  is  the  fact  that  they  converge  to  a  Gaussian  as 
their  order  tends  to  infinity.  Unser  et  oil.  [16]  derived  the  Gaussian  approximation 

Pp(x)  -  -75=  e 

yzircjp 

where  ap  =  Ratio  ,  for  example,  is  already  well  below  0.2%. 

Denoting  by  Sp  a  spline  function  space  spanned  by  {Pp(x  -  m)}  for  p  odd  and  by 
{Pp(x  m)}  for  p  even  (subscript  in  S\  refers  to  the  fact  that  spline  functions  have 
knots  at  integers),  spline  spaces  form  a  nested  sequence 

. . .  C  Sp2i  C  ...  C  Sp  C  Sp  C  Sv2_  j  C  ...  C  Sp2-i  C - By  orthogonalizing  this  basis 

functions  a  multiresolution  analysis  of  IP  [R) .  from  which  the  Battle-Lemarie  wavelet  bases 
stem,  can  be  built  [17].  Here,  we  will  not  pursue  constructions  of  orthogonal, 
semi-orthogonal,  or  biorthogonal  spline- wavelets  any  further— reader  looking  for  a  detailed 
treatment  of  this  subject  may  find  a  good  starting  point  in  books  [18,  19]. 

2.1.2  B-Spline  Signal  Interpolations 

Unser  et  al.  developed  a  fast  digital  filtering  scheme  for  B-spline  signal  processing  [20]. 
They  defined  a  discrete  B-spline  of  order  p  and  expansion  factor  (spacing  between  knots)  m 
as 

Pp  ’  n,  m  £  Z. 

Henceforth,  if  the  distance  between  knots  equals  one,  we  will  write  bp{n)  instead  of  bPti(n). 
Interpolation  of  a  discrete  signal  s(n)  6  12{Z)  by  sp(x)  6  Sp  using  central  B-splines 

OO 

s(n)  =  sp(x )  =  c(i)Pp(x  ~  i)  ,  (1) 

x—n  *=— oo  x=n 

can  now  be  written  as  a  convolution  sum 

s(n)  =  c*bp(n).  (2) 

If  s(n)  are  samples  of  a  function  s(x)  bandlimited  to  [-7r,7r]  (i.e.,  the  support  of  its  Fourier 
transform  S(u)  is  in  [— 7r,  7r]),  it  can  be  shown  that  sp(a;)  ->■  s(a:)  as  p  ->  oo  [21,  22]. 

In  [20],  they  refer  to  a  linear  operator  by  which  B-spline  coefficients  c{n)  can  be  obtained 
from  samples  s(n)  as  a  “direct  B-spline  transform.”  Equation  (2)  therefore  represents 
“indirect  B-spline  transform”  of  a  sequence  (c(n)}. 

After  taking  the  ^-transform  of  (2),  the  direct  B-spline  filters  are  found  to  be 

bp 1  (n)  =  1{[-Bp(^)]  1}.  Since  B-spline  functions  Pp{x)  are  compactly  supported,  indirect 
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Table 


P 


Transfer  functions  of  direct  B-spline  filters  for  orders  from  0  to  9. 
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z+6+z- 1 
6 

z+4+z-1 

384 


z2+76z+230+76z-1+z~2 

_ 120 _ 

22-f26z+66+26z“1+z“2 

46080 


^3+722z2+10543>2+23548+10543^-1+7222-2+2-3 

_ 5040 _ 

23+120z2+U9l2+2416+119l2-1+120;z-2+2-3 

_ 10321920 _ 

*4+6552z3+331612z2+2485288z+4675014+2485288;z-1+331612z-2+6552z“3+;z-4 

_ 362880 _ 

_ £4+50223+14608z2+88234z4-156190+88234z-1-Tl4608z-2-l-502;z--3+z-4 


B-spline  filters  bp(n)  are  finite  impulse  response  (FIR)  filters,  while  direct  B-spline  filters 
b~l{n)  are  infinite  impulse  response  (HR)  filters.  Aldroubi  et  al.  [21]  showed  that  HR 
filters  bp  1  (n)  are  stable  (i.e.,  the  region  of  convergence  of  B~l(z)  includes  the  unit  circle 
[23])  for  any  order  p.  Note  that  both  indirect  and  direct  B-spline  filters  are  symmetric, 
which  follows  from  the  fact  that  central  B-splines  /3p(x)  are  symmetric. 

Table  1  shows  the  ^-transforms  of  direct  B-spline  filters  for  the  first  ten  orders.  We 
postpone  the  discussion  on  implementation  details  of  B-spline  filters  until  Sections  2.2.8 
and  2.2.9. 

Instead  of  using  B-spline  interpolation  as  given  by  (1)  it  is  sometimes  convenient  to  express 
the  interpolating  function  sp(x)  in  terms  of  discrete  samples  s(n) 

OO 

sp(x)  =  £  s{i)r)p(x-i),  (3) 

i~ — OO 

where  r}p(x)  is  the  cardinal  spline  of  order  p.  In  the  frequency  domain,  cardinal  splines 
converge  to  an  ideal  lowpass  filter  with  cutoff  frequency  tt  (i.e.,  r)p(x)  ->•  sine  (a;))  as  p  tends 
to  infinity  [21,  24],  which  establishes  the  asymptotic  equivalence  with  Shannon’s  sampling 
theorem  [25]. 
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(4) 


Using  (1)  and  (2)  with  (3)  cardinal  splines  can  be  related  to  B-splines: 

OO 

V pix)  bp  —  *)• 

i=— oo 

Cardinal  splines  r]p(x)  and  rjp(u)  for  p  G  {0, 1, 2, 3, 4}  are  shown  in  Figure  2. 

2.1.3  B-Spline  Signal  Approximations 

Central  B-splines  are  also  simple  to  use  when  the  goal  is  function  approximation. 
Least-squares  B-spline  approximation  of  s(z)  G  L2(R)  is  achieved  by  computing  the 
orthogonal  projection  of  this  function  onto  Sp.  We  have 


OO 


*(*)  =  sp(x)  =  J2  d(i)(3p(x  -  i) 

i~—oo 

(5) 

with 

d(i)  =  (s(x)Jp(x-i)), 

(6) 

where 

O  00 

(3p(x)  =  ^2  b2p+i(i) f3p(x  —  i) 

i=- oo 


is  the  dual  spline  of  order  p  [24].  Spline  functions  (5p(x)  and  f3p(x)  satisfy  the 
biorthogonality  condition 

o 

(Pp(x  -  m),/3p(x  -  n))  =  5(m  -  n),  m,neZ. 

O 

(Note  that,  since  both  Pp(x)  and  f3p(x)  form  a  basis  of  Sp,  they  can  be  interchanged  in  (5) 
and  (6).) 

o  o 

Figure  3  shows  functions  f)p{x)  and  their  Fourier  transforms  /3p(w)  for  p  G  {0, 1,  2, 3, 4}. 

An  interesting  alternative  to  the  minimum  L2-norm  (i.e.,  least-squares)  approximation  of  a 
signal  is  obtained  by  computing  oblique  instead  of  orthogonal  projection  of  the  signal  onto 
the  spline  function  space.  Unser  and  Aldroubi  proposed  an  independent  specification  of  the 
sampling  and  approximation  spaces  [26]:  a  linear  operator  maps  coefficients  of  the  input 
signal  expansion  over  sampling  space  basis  into  the  coefficients  of  the  approximation  space 
basis  expansion  from  which  the  signal’s  projection  onto  the  approximation  space  is 
recovered.  Constraining  the  entire  system  to  be  linear,  shift-invariant  for  integer 
translations,  and  consistent  (i.e.,  the  system  acts  as  an  identity  operator  for  functions  that 
belong  to  the  approximation  space),  the  obtained  solution  for  the  signal  approximation  is 
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the  projection  of  a  signal  onto  the  approximation  space  perpendicular  to  the  sampling 
space  [26].  Analogously  to  (5)  and  (6)  this  projection  can  be  expressed  as 

00 

»r(®)  =  ^  a(i)/3r(x  -  i)  (7) 

i~~oc 

with 

=  Qn  *  (s{x),Ps(x  ~  *)>» 

where  is  the  convolution  inverse  of  the  cross-correlation  sequence 

912(0  =  (P»(x  -  0.  Pr(x))  =  bs+r+i(i). 

When  the  sampling  space  Ss  and  the  approximation  space  Sr  are  identical  (i.e.,  s  =  r),  an 
orthogonal  projection  given  by  (5)  and  (6)  results.  (Note  that  the  described  oblique 
projection  is  not  restricted  to  spline  function  spaces — the  only  requirement  is  that  both 
sampling  space  basis  and  approximation  space  basis  are  Riesz  bases  of  the  corresponding 
function  spaces  [26].) 

Signal  approximation  (7)  is  particularly  attractive  in  situations  where  the  sampling  space  is 
given  a  priori  (e.g.,  by  the  impulse  response  of  the  acquisition  device  [26])  or  when  such  an 
approximation  is  close  to  the  optimal  least-squares  solution  but  simpler  to  implement  than 
the  orthogonal  projection  (e.g.,  [27]). 

2.2  Steerable  Dyadic  Transforms 

In  this  section,  we  derive  transforms  that  will  enable  efficient  directional  processing  of 
mammograms.  Section  2.2.1  deals  with  problems  due  to  the  lack  of  translation  and 
rotation  invariance  of  orthogonal  and  biorthogonal  wavelet  transforms.  In  Section  2.2.2, 
the  one-dimensional  discrete  dyadic  wavelet  transform  [11]  is  augmented  to  obtain  a  strong 
foundation  for  derivations  of  two-dimensional  transforms.  Connections  to  reconstruction 
from  edges,  and  to  scale-space  filtering  are  mentioned  in  Section  2.2.3.  Section  2.2.4 
presents  a  direct  extension  of  the  one-dimensional  transform  to  two  dimensions.1  The 
concept  of  steerability  is  explained  in  Section  2.2.5,  and  then  used  for  derivations  of  a 
steerable  dyadic  wavelet  transform  and  multiscale  spline  derivative-based  transform  in 
Sections  2.2.6  and  2.2.7,  respectively.  Sections  2.2.8  and  2.2.9  describe  how  to  exploit 
symmetries  of  signals  and  filters  for  a  fast  implementation  of  the  transforms. 

1For  extensions  to  higher  dimensions,  please  refer  to  [12,  28]. 
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2.2.1  Shortcomings  of  Traditional  Methods  of  Wavelet  Analysis 

Analyzing  images  across  multiple  scales  and  resolution  has  become  a  powerful  tool  for 
solving  compelling  problems  in  computational  vision,  image  processing,  and  pattern 
recognition.  Wavelet  theory  encompasses  multiscale  and  multiresolution  representations, 
such  as  subband  filtering  [29],  image  pyramids  [30],  and  scale  space  filtering  [31],  into  a 
unified  mathematical  framework.  In  the  area  of  image  processing,  there  remain  few 
research  areas  to  which  wavelet  analysis  has  not  been  applied.  For  example,  problems  in 
image  compression,  denoising,  restoration,  enhancement,  registration,  fusion,  segmentation, 
and  analysis,  have  all  been  approached  with  distinct  kinds  of  wavelet  processing. 

Though  ubiquitous,  wavelet  analysis  is  not  without  problems  of  its  own.  Lack  of 
translation  invariance,  one  of  the  major  problems  of  the  wavelet  transform  [17],  is  in 
multiple  dimensions  accompanied  with  lack  of  rotation  invariance. 

Wavelet  transform  in  its  most  commonly  used  orthogonal  or  biorthogonal  forms  is  not 
translation  and  rotation-invariant.  By  translation-invariant  transform,  we  mean  a 
transform  that  commutes  with  a  translation  operator.  Since  we  will  deal  primarily  with 
discrete  transforms  in  this  work,  we  constrain  the  translation  parameter  to  integer 
multiples  of  a  sampling  period. 

Lack  of  translation  invariance  of  the  discrete  wavelet  transform  is  illustrated  in  Figure  4. 
Here,  we  can  clearly  see  how  a  translation  of  the  input  signal  by  one  sample  results  in  a 
completely  different  set  of  transform  coefficients  (orthogonal  wavelet  DAUB42  [17]  was  used 
in  this  experiment). 

Noninvariance  under  translations  of  an  orthogonal  and  biorthogonal  wavelet  transform  is 
due  to  lower  sampling  density  at  coarser  scales.3  A  straightforward  way  of  dealing  with  this 
problem  is  to  construct  a  redundant  transform  by  using  the  same  sampling  frequency  for 
the  input  signal  and  all  scales  of  the  transform.  A  filter  bank  implementation  of  such  a 
transform,  called  “algorithme  a  trous”  [32],  is  based  upon  the  fact  that  downsampling 
followed  by  filtering  is  equivalent  to  filtering  with  the  upsampled  filter  before  the 
downsampling,  as  shown  in  Figure  5. 

Lack  of  rotation  invariance  is  another  shortcoming  of  traditional  (i.e.,  orthogonal  and 
biorthogonal)  wavelet  techniques.  In  defining  rotation  invariance,  we  are  a  bit  less  strict 
than  with  translation  invariance.  We  do  not  require  that  the  transform  commutes  with  a 
rotation  operator  here.  Even  in  the  case  of  a  simple  filtering,  this  would  limit  us  to 

2The  number  in  DAUB4  refers  to  twice  the  order  of  the  wavelet  (i.e.,  two  in  this  case). 

3In  practice,  since  analysis  is  performed  over  a  finite  range  of  scales,  a  discrete  wavelet  transform  is 
translation-invariant  by  translations  determined  by  the  coarsest  scale  (e.g.,  sixteen  samples  for  the  analysis 
from  Figure  4)  [17]. 
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Figure  5:  Filter  bank  implementation  for  (a)  a  discrete  wavelet  transform  and  (b)  “algorithme 
a  trous”  decompositions  for  three  levels  of  analysis. 

circularly  symmetric  filters  only.  Our  requirement  for  analysis  is  a  transform  that  enables 
rotation-invariant  processing.  As  an  example  of  such  a  transform,  let  us  consider  filtering 
with  the  first  derivative  of  a  two-dimensional  Gaussian  probability  density  function  in  two 
directions,  specifically,  along  x  and  along  y- axis.  By  linearly  combining  the  results  of 
filtering  in  these  two  directions,  filtering  with  the  first  derivative  of  a  Gaussian  in  any 
direction  can  be  computed.  This  fact  was  used  by  Canny  [33]  for  edge  detection.  A 
determined  edge  direction  rotates  as  an  input  image  is  rotated. 

After  choosing  the  fundamental  properties  of  the  transform,  one  must  decide  upon  the 
basis  functions  to  be  applied.  For  our  studies,  we  selected  basis  functions  that  well 
approximated  derivatives  of  a  Gaussian,  because  (1)  the  Gaussian  probability  density 
function  is  optimally  concentrated  in  both  time  and  frequency  domain,  and  thus  suitable 
for  time-frequency  analysis,  (2)  higher  order  derivatives  of  a  Gaussian  can  be,  similar  to 
the  first  derivative,  used  for  rotation-invariant  processing  [34],  and  (3)  the  Gaussian 
function  generates  a  causal  (in  a  sense  that  a  coarse  scale  depends  exclusively  on  the 
previous  finer  scale)  scale  space  [35].  The  last  property  makes  possible  scale-space 
“tracking”  of  emergent  features. 

2.2.2  1-D  Discrete  Dyadic  Wavelet  Transform  Revisited 

A  discrete  wavelet  transform  is  obtained  from  a  continuous  representation  by  discretizing 
dilation  and  translation  parameters  such  that  the  resulting  set  of  wavelets  constitutes  a 
frame.  The  dilation  parameter  is  typically  discretized  by  an  exponential  sampling  with  a 
fixed  dilation  step  and  the  translation  parameter  by  integer  multiples  of  a  fixed  step  [17]. 
Unfortunately,  the  resulting  transform  is  variant  under  translations,  a  property  which 
makes  it  less  attractive  for  image  analysis. 
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As  we  have  already  mentioned  in  Section  2.2.1,  sampling  the  translation  parameter  with 
the  same  sampling  period  as  the  input  function  to  the  transform  results  in  a 
translation-invariant,  but  redundant  representation.  The  dyadic  wavelet  transform 
proposed  by  Mallat  and  Zhong  [11]  is  one  such  representation.  Let  us  begin  with  a  brief 
review  of  properties  of  the  dyadic  wavelet  transform  as  described  in  [11],  but  included  here 
for  completeness. 

The  dyadic  wavelet  transform  of  a  function  s(x)  e  L2(R)  is  defined  as  a  sequence  of 
functions 

{bLms(a;)}me^,  (8) 

where  Wms(x )  =  s  *  ij)m{x),  and  ipm(x)  =  2~mip(2~mx)  is  a  wavelet  ip(x )  expanded  by  a 
dilation  parameter  (or  scale)  2m.  Note  the  use  of  convolution  instead  of  an  inner  product. 
To  ensure  coverage  of  the  frequency  axis  the  requirement  on  the  Fourier  transform  of 
ipm(x)  is  the  existence  of  A\  >  0  and  B i  <  oo  such  that 

OO 

Ai  <  £  \xp(2muj)\2  <  Bl 

171— —  OO 

is  satisfied  almost  everywhere.  The  constraint  on  the  Fourier  transform  of  the  (nonunique) 
reconstructing  function  x(x)  is 

OO 

E  i’(2mu)x(2muJ)  =  l. 

771— —  OO 

A  function  s(x)  can  then  be  completely  reconstructed  from  its  dyadic  wavelet  transform 
using  the  identity 

OO 

s(a;)  =  ^  )  iFTOs  *  Xm(x)t 

m=— oo 

where  Xm(x)  =  2~mx{2~mx). 

In  numerical  applications,  processing  is  performed  on  discrete  rather  than  continuous 
functions.  When  the  function  to  be  transformed  is  in  the  discrete  form,  the  scale  2m  can  no 
longer  vary  over  all  m  £  Z.  Finite  sampling  rate  prohibits  the  scale  from  being  arbitrarily 
small,  while  computational  resources  restrict  the  use  of  an  arbitrarily  large  scale.  Let  the 
finest  scale  be  normalized  to  1  and  the  coarsest  scale  set  to  2M. 

The  smoothing  of  a  function  s(x)  €  L2(R)  is  defined  as 

Sms(x)  =  S  *  4>m{x)i 

where  (fim(x )  =  2~m<j){2~mx )  with  m  £  Z,  and  (f>{x)  is  a  smoothing  function  (i.e.,  its 
integral  is  equal  to  1  and  <f>(x)  0  as  |a;|  ->■  oo). 
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(9) 


In  [11],  a  real  smoothing  function  <f>(x)  was  selected,  whose  Fourier  transform  satisfied 

OO 

I^MI2=  EV’(2"V)x(2mw). 

m~\ 

In  addition,  it  was  shown  that  any  discrete  function  of  finite  energy  s(n)  €  l2(Z)  can  be 
written  as  the  uniform  sampling  of  some  function  smoothed  at  scale  1,  i.e.,  s(n)  =  Sof(n), 
where  f(x )  G  L2(R )  is  not  unique.  Thus,  the  discrete  dyadic  wavelet  transform  of  s(n)  for 
any  coarse  scale  2M  was  defined  as  a  sequence  of  discrete  functions 


{ S M f  {Tl  +  s),  {Wmf  (n  +  ^)}m£[l,M]}ng^i 
where  s  is  a  ip(x)  dependent  sampling  shift. 

The  above  initialization  s(n )  =  SQf(n)  is  rather  standard  in  the  discrete  wavelet  transform 
computation  [17],  although  it  yields  correct  results  (i.e.,  the  discrete  wavelet  transform  is 
equal  to  the  samples  of  its  continuous  counterpart)  only  when  s(n)  =  S'os(n).  Here,  we  will 
concentrate  on  wavelets  which  are  derivatives  of  spline  functions  and  this  will  lead  us  to  a 
simple  initialization  procedure  [36]  that  alleviates  the  above  problem. 

For  a  certain  choice  of  wavelets  the  discrete  dyadic  wavelet  transform  can  be  implemented 
within  a  fast  hierarchical  digital  filtering  scheme.  Next,  we  shall  summarize  the  relations 
between  filters,  wavelets,  and  smoothing  functions. 

First,  let  us  introduce  a  real  smoothing  function  <p(x)  such  that  (9)  can  be  rewritten  as4 

OO 

${u)  <p{uj)  =  ’4)(2rnuj)  x(2mw),  (10) 

m=0 

and  let  us  set  4>{x)  =  (5P (x)  (i.e.,  we  restrict  ourselves  to  wavelets  which  are  spline 
functions). 

Computing  (10)  for  the  finest  two  scales  shows  that 


xM  =  AH  -  A(H  <£(2 oo). 

/3p(2u)  can  be  related  to  (3p(uj)  by  expressing  /3p(2a>)  as  (cf.  Proposition  1  of  [36]) 

P+1  7 sin  (f)X  P+1 


(11) 


Pp(  2w) 


2p+i 


sin(u;) 

sin  (l). 


U 

2 


and  using  E"=0  ej(mtJ+e)  =  s'n^^)  ^  ej^+e): 


a»(2 w) = (cos  AH- 


(12) 


4Note  that  the  sum  index  determines  the  range  of  scales  of  the  discrete  transform:  using  (9)  we  have 
^(2w)  and  x(2cj)  at  the  finest  scale  of  the  transform,  while  for  (10)  we  get  $(u)  and  x(u)- 
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(Note  that  a  relation  similar  to  (12)  can  be  derived  for  integer  scales  provided  that  the 
dilation  parameter  and  the  order  p  are  not  both  even  [36].) 

Let  F(u)  be  a  digital  filter  frequency  response  and  let  Fs( u)  —  e^wsF(u). 

If  we  choose 


=  G-a(u)pp(u), 

(13) 

<P(  2w)  =  La(u)ip((o), 

(14) 

XW)  =  Ks(u)(p(u}), 

(15) 

HH  =  (cos  (|))’’+1 , 

(16) 

where  s  €  {0,  |}  is  a  filter  dependent  sampling  shift  needed  for  g(n),  l(n),  k(n),  and  h(n ) 
to  be  FIR  filters,  and  insert  Equations  (12)— (16)  into  (11),  we  observe  the  relation  between 
frequency  responses  of  the  filters 

G(u)K(u>)  +  H(oo)L(u>)  =  1.  (17) 

Similar  to  orthogonal  and  biorthogonal  discrete  wavelet  transforms,  the  discrete  dyadic 
wavelet  transform  can  be  implemented  within  a  hierarchical  filtering  scheme.  To  derive 
such  a  digital  filtering  scheme,  let  us  assume  that  s(w)  from  (8)  is  bandlimited  to  [-7r,7r]. 
Using  Shannon’s  sampling  theorem  [25]  and  (13)  in  the  definition  of  the  dyadic  wavelet 
transform  (8)  with  m  =  0  shows 

/OO  °°  00 

s(z)sinc(t  —  i)  ]T  g_a(m)(3p(x  -  t  -  m)  dt. 

"°°  —oo  oo 

By  making  use  of  the  fact  that  the  cardinal  spline  functions  tend  to  the  sine  function  as 
their  order  r  approaches  infinity,  and  employing  (4)  we  can  write 

feM  *  S{u)B;\u)pr{u)PP{u)G-s{u), 
or,  by  using  (12)  and  (16), 

m— 1 

F{Wms(x)\x=n}  ~  S(u>)  B;\uj)  Bp+r+i(u>)  G_s(  2"*u;)  ft  H-tfu,).  (18) 

i= 0 

Equation  (18)  entirely  specifies  the  discrete  dyadic  wavelet  transform  decomposition,  while 
the  reconstruction  follows  from  (11)— (16).  Three  levels  of  a  filter  bank  implementation  are 
shown  in  Figure  6.  (Note  that  the  initialization  is  the  same  as  the  one  proposed  in  [36]  and 
that  except  for  the  prefiltering  and  postfiltering,  this  scheme  is  implementing  “algorithme  a 
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I  I  GJ«> 


^((O) 

Figure  6:  Filter  bank  implementation  of  a  one-dimensional  discrete  dyadic  wavelet  transform 
decomposition  (left)  and  reconstruction  (right)  for  three  levels  of  analysis. 

trous.”)  Noninteger  shifts  at  scale  1  for  filters  with  s=|  are  rounded  to  the  nearest  integer. 


We  will  now  perform  a  simple  experiment  which  will  illustrate  the  difference  between  the 
implementation  of  the  discrete  dyadic  wavelet  transform  as  originally  proposed  in  [11]  (i.e., 
without  prefiltering  and  postfiltering)  and  the  one  from  Figure  6. 

Let  s(x)  =  sinc(a;),  p  =  2,  and  g(n)  =  2 S(n  +  1)  —  2 S(n)  (this  particular  choice  for  p  and 
g(n)  results  in  the  same  wavelet  as  was  used  by  Mallat  and  collaborators  [37,  11]).  The 
dyadic  wavelet  transform  of  s(:r)  at  a  scale  2m  (8)  in  the  frequency  domain  is  then 

wZs{u)  =  G_s( 2mu)  &(2mu;)  rect^)  .  (19) 

The  Fourier  transform  of  the  discrete  dyadic  wavelet  transform  of  s(n)  =  S(n)  at  a  scale  2m 
using  spline  based  initialization  follows  from  (18) 


F{Wms(n)}  =  B;1(u)Br+3 


m- 1 


(w)G-.eru) 

i=0 


(20) 


while  the  one  using  the  algorithm  from  [11]  equals 


f{Wms(n)} 


m— 1 


G-s( 2mu)  JJ  H-S(2iu>). 


i=0 


(21) 


In  Figure  7  a  comparison  of  magnitudes  of  (20)  and  (21)  versus  (19)  is  shown:  in  Figure 
7(a)  magnitudes  of  (19)  (solid)  and  (21)  (dashed)  are  plotted  for  m  e  {0, 1, 2,3},  while  the 
dashed  curves  in  7(b)  represent  magnitudes  of  (20)  with  r= 5. 

By  choosing  the  appropriate  order  r,  (20)  can  approximate  (19)  in  the  interval  [-7 r,  7r] 
arbitrarily  good,  while  originally  proposed  (21)  has  troubles  at  finer  scales.  Mallat  and 
Zhong  [11]  noticed  that  there  was  a  problem  with  their  discrete  transform  computation, 
and  introduced  a  set  of  constants  associated  with  the  discrete  transform  coefficients  at 
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0  x  -T  0 

(a)  (b) 


Figure  7:  (a)  Fourier  transform  magnitudes  of  the  dyadic  wavelet  transform  of  s(x )  =  sinc(a;) 
(solid)  and  the  originally  proposed  discrete  dyadic  wavelet  transform  [11]  of  s(n)  =  5(n) 
(dashed),  (b)  Fourier  transform  magnitudes  of  the  dyadic  wavelet  transform  of  s(x )  (solid) 
and  the  discrete  dyadic  wavelet  transform  using  quintic  splines  for  interpolation  of  s(n) 
(dashed) . 
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dyadic  scales.  They  chose  the  values  of  constants  such  that  the  transform  coefficient 
modulus  maxima  remained  constant  over  all  dyadic  scales  for  a  step  edge  input  signal.  In 
relation  to  Figure  7(a)  this  is  equivalent  to  multiplying  F{Wms{n)}  by  a  distinct  constant 
for  each  m.  Clearly,  this  can  improve  over  the  situation  depicted  by  Figure  7(a),  but  can 
still  not  compete  with  the  described  spline  based  initialization. 

Next,  we  will  choose  filters  in  the  filter  bank  implementation  of  the  discrete  dyadic  wavelet 
transform.  As  already  mentioned,  we  are  interested  in  wavelets  which  are  derivatives  of 
spline  functions.  G(u ;)  in  (13)  is  therefore  the  Fourier  transform  of  the  difference  operator 
centered  around  —  s  (cf.  Property  4  in  Section  2.1): 

G(w)  =  e*"(2isin(|))',  (22) 


where  d  is  the  order  of  the  derivative  and  the  sampling  shift  for  this  filter  is  s  =  dm°d2. 
Since  H(ui)  was  already  given  by  (16),  the  remaining  two  filters  to  be  determined  are  L(cu) 
and  K(u).  Both  of  them  are  (as  is  true  for  ip(x)  and  x(z))  nonunique. 

If  we  choose  L(uj)  such  that  we  can  express  K(uj)  in  terms  of  a  finite  geometric  series 
having  the  smallest  number  of  elements  for  an  arbitrary  p,  we  get 


and 


l^J 

LH  =  £  (-1) 


orst-en" 


(23) 


(24) 


where  [xj  denotes  the  largest  integer  smaller  than  x,  the  sampling  shift  for  L(u)  is  the 
same  as  the  one  for  H(u)  (i.e.,  s=  (?>fl)^nod2))  ancj  the  sampling  shift  for  K(u)  is  the  same 
as  the  one  for  G(lo). 

Note  that  Equations  (16)  and  (22)-(24)  work  fine  from  the  mathematical  point  of  view, 
but  in  practice  the  reconstruction  may  become  cumbersome  when  both  p  and  d  are  large 
(the  lengths  of  impulse  responses  h(n),  g(n),  l(n),  and  k(n)  are  p+ 2,  d+ 1, 

(p+l)(d—  (d+l)mod2)+l,  and  pd+(p+l)(dmod2)+l,  respectively,  while  for  the  frequency 
responses  of  the  decomposition  filters  we  observe  that  limp-**,  |/f_s(u>)|  =  5u(u  +  2nn)  and 
limd_>O0(2j)_d  |G_s(o;)|  =  8u(oo  +  (2n+l)7r)  with  n  6  Z ). 

It  is  also  worth  noting  that  K(u>)  is  a  lowpass  filter  when  p  is  even  (i.e.,  the  reconstruction 
function  x(^)  is  a  wavelet  only  for  p  odd). 

Tables  2,  3,  and  4  list  impulse  responses  of  the  four  filters  for  p  G  {0, 1, 2}  and  d  G  {1, 2, 3}, 
while  Figure  8  shows  wavelets  ip(x)  =  -  for  the  same  values  of  p  and  d.  Wavelets 
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Table  2:  Impulse  responses  /i(n),  g(n),  /(re),  and  fc(n)  for  p= 0  and  d  G  {1,2,3}. 


d=  1 

d  =  2 

d  =  3 

n 

h(n ) 

s(«) 

Kn) 

k(n) 

9(n ) 

/(n) 

k(ri) 

P(«) 

/(n) 

k(n) 

-2 

1 

-1 

0.5 

1 

1 

-3 

-0.125 

0 

0.5 

-1 

0.5 

-0.25 

-2 

0.5 

-0.25 

3 

0.625 

0.0625 

1 

0.5 

0.25 

1 

0.5 

-1 

0.625 

-0.0625 

2 

-0.125 

Table  3:  Impulse  responses  h 


n),  g(n),  Z(n),  and  A:(n)  for  p=l  and  d  e  {1, 2, 3}. 


h[n) 

d 

=  1 

d 

=  2 

n 

l(n) 

9(n) 

k(n) 

P(«) 

k(n) 

-i 

1 

Bajiiwa 

0 

0.5 

-1 

-0.375 

1 

2 

0.25 

0.3125 

0.0625 

-0.0625 

CO 

II 

n 

h(n ) 

g(n) 

l(n) 

k(n) 

-0.015625 

i 

-0.09375 

0.00390625 

0.25 

-3 

0.265625 

0.04296875 

0 

0.5 

3 

0.6875 

0.1015625 

1 

0.25 

-1 

0.265625 

-0.1015625 

2 

-0.09375 

-0.04296875 

3 

-0.015625 

-0.00390625 
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(b) 


(c) 


Figure  8:  (a)  Wavelets  (a)  i^(x)  =  d^p^fx\  (b)  ip(x)  =  >  and  (c)  ip(x)  =  for 

p=0  (dashdotted),  p= 1  (solid),  and  p= 2  (dashed). 
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Table  4:  Impulse  responses  /t(n),  g(re),  l(n),  and  kjn )  forp=2  and  d  e  {1, 2, 3} 


d  = 

1 

2 

n 

h(n) 

HQl 

l{n) 

k(n) 

9(n) 

l(n) 

k(n) 

-2 

-1 

0.125 

0.375 

■ 

0.125 

-0.015625 

-0.109375 

1 

0.125 

-0.015625 

-0.125 

0 

0.375 

0.375 

-0.34375 

-2 

0.375 

-0.46875 

1 

0.125 

0.375 

0.34375 

1 

0.375 

-0.125 

2 

3 

0.125 

0.109375 

0.015625 

0.125 

-0.015625 

d  =  3 

n 

h(n) 

9{n) 

l{n) 

k(n) 

-4 

-0.001953125 

0.000244140625 

-3 

-0.017578125 

0.003662109375 

-2 

B 

1 

-0.0703125 

0.0263671875 

-1 

s 

-3 

0.085937 

0.0908203125 

0 

eM 

3 

0.50390625 

0.13037109375 

1 

0.125 

-1 

0.50390625 

-0.13037109375 

2 

0.0859375 

-0.0908203125 

3 

-0.0703125 

-0.0263671875 

4 

-0.017578125 

-0.003662109375 

5 

-0.001953125 

-0.000244140625 

from  this  family  have  a  support  of  length  d+p+1,  regularity  order  p,  and  are  either 
symmetric  or  antisymmetric. 

As  already  discussed,  wavelets  with  p—2  and  d= 1  from  a  family  of  wavelets  with  p  even 
and  d=l  were  used  in  [37,  11],  whereas  filters  with  p= 1  and  d= 2  from  a  family  of  filters 
with  p  odd  and  d=2  were  employed  by  Laine  and  collaborators  [10,  8,  7].  Here  described 
transform  puts  no  restrictions  on  the  choice  of  p  or  d  whatsoever,  and  uses  a  better 
initialization  procedure  than  the  one  originally  proposed  in  [11]. 

2.2.3  Remarks 

1.  As  mentioned  in  Section  2.2.1,  the  translation  invariance  property  of  the  discrete  dyadic 
wavelet  transform  is  due  to  the  fact  that  the  translation  parameter  is  sampled  with  the 
same  sampling  period  as  the  input  signal  over  all  scales.  When  comparing  Figure  6  to 
traditional  filter  bank  implementations  of  an  orthogonal  and  biorthogonal  wavelet 
transform,  we  observe  that  each  subband  shown  in  Figure  6  is  retained  at  its  original 
density  rather  than  being  critically  sampled.  Note,  that  the  filters  at  distinct  scales  are  not 
the  same  as  in  orthogonal  or  biorthogonal  cases  of  analysis  [17]. 

The  discrete  dyadic  wavelet  transform  is  highly  redundant.  To  obtain  a  more  parsimonious 
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representation  we  may  sample  the  translation  parameter  in  a  translation-invariant  manner 
(e.g.,  sampling  a  function  at  its  extrema  or  extrema  of  its  derivatives).  Of  importance  is, 
whether  or  not  it  remains  possible  to  reconstruct  the  original  signal  from  this  subset  of 
transform  coefficients. 

In  [37,  11]  a  reconstruction  algorithm  is  presented  that  reconstructs  an  approximation  of  a 
signal  given  the  position  of  local  maxima  of  wavelet  coefficients  modulus  and  the  values  of 
wavelet  coefficients  at  each  corresponding  location.  Wavelets  with  d— 1  were  used  so  that 
the  locations  of  maxima  corresponded  to  inflection  points  of  the  original  signal  smoothed 
at  dyadic  scales. 

2.  For  cases,  when  the  wavelet  ip(x)  from  Section  2.2.2  is  equal  to  either  the  first  (d=  1)  or 
the  second  (d-  2)  derivative  of  a  smoothing  function  d(x)  =  /3p+d(x)  the  dyadic  wavelet 
transform  of  a  function  s(x )  €  L2(R )  can  be  written  as 

Wms(x )  =  s  *  ^2md  ~  j  (a;)  =  2mrf-j^j(s  *  0m)(x)  meZ,  de{l,2}. 

Depending  on  the  wavelet  selected  for  analysis,  the  wavelet  transform  Wms(x)  is 
proportional  either  to  the  first  or  second  derivative  of  s(x)  smoothed  by  9m.  By  recording 
the  zero-crossings  of  Wms(x)  for  d  =  2,  a  scale-space  image  of  the  representation  similar  to 
[31]  can  be  obtained.  The  only  differences  are  that  curves  in  the  scale-space  plane  are 
computed  for  dyadic  scales  only  and  that  dm  is  a  close  approximation  of  a  Gaussian 
instead  of  a  true  Gaussian. 

2.2.4  2-D  Discrete  Dyadic  Wavelet  Transform  Revisited 

The  dyadic  wavelet  transform  of  a  function  s(x,y )  e  L2(R2)  is  defined  as  a  set  of  functions 
[11] 

{W^s(x,y),Wls(x,y)}meZ,  (25) 

where  W^s(x,  y)  =  s  *  ip^x,  y)  for  i  =  {1, 2}  and  y)  =  2-2mipi(2~mx,  2 ~my)  are 
wavelets  ipl(x,y)  expanded  by  a  dilation  parameter  2m. 

To  ensure  coverage  of  the  frequency  space  there  must  exist  an  A2  >  0  and  B2  <  oo  such 
that 

oo  2 

^2  <  £  £1  i’i(2muJx,2muy)\2  <  b2 

m=— oo  {— 1 

is  satisfied  almost  everywhere.  If  (nonunique)  functions  yx(a;,y),  y2(a :,y)  are  chosen  such 
that  their  Fourier  transforms  satisfy 

oo  2 

E  E  2”w„)  Xi(2mwI,  2’”w„)  =  1, 
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the  function  s(x,y)  may  be  reconstructed  from  its  dyadic  wavelet  transform  by 

00  2 

«(*,»)=  £  £<***!»(*.»). 

m=— oo  i—1 

where  xLix,v)  =  2-2mxi{2~mx,2 ~my). 

However,  when  processing  discrete  functions  the  scale  2m  may  no  longer  vary  over  all 
m  G  Z.  Let  the  finest  scale  be  normalized  to  1  and  the  coarsest  scale  set  to  be  2M.  Let  us, 
similar  to  [11],  introduce  a  real  smoothing  function  (f)(x,y )  such  that  its  Fourier  transform 
satisfies 

oo  2 

l*(w„w,)|a  =  £  ^  2muy)  y?(2mux,  2mwy).  (26) 

m=0  2—1 

Here,  as  in  one  dimension,  a  finite  energy  discrete  function  ( s(nx,ny )  G  l2(Z2))  can  be 
written  as  the  uniform  sampling  of  some  function  smoothed  at  scale  1:  s(nx,  ny) 

—  Sof(nx,ny),  where  f(x,y)  €  L2(R2)  is  not  unique,  and  Smf(x,  y)  =  f*  (f>m(x,  y).  This 
led  Mallat  and  Zhong  [11]  to  a  two-dimensional  analog  of  the  one-dimensional  definition  of 
the  discrete  dyadic  wavelet  transform:5 

{SM-if(nx  +  s,ny  +  s ),  {W^f(nx  +  s,  ny  +  s),  W^f{nx  +  s,ny  + 

We  will  use,  as  in  Section  2.2.2,  a  spline-based  initialization  procedure. 

To  implement  a  multidimensional  discrete  dyadic  wavelet  transform  within  a  fast 
hierarchical  digital  filtering  scheme,  the  wavelets  were  chosen  to  be  separable  products  of 
one-dimensional  functions: 


^(x,y)  =  ^(x)  <f>(y),  (27) 

^2(x,y)  =  0(y)  <j>{x),  (28) 

where  <j>(x)  and  ^( x )  were  chosen  as  described  in  Section  2.2.2  (i.e.,  4>(x)  =  f3p(x)  and  0(u;) 
specified  by  (13)). 

From  (27),  (28),  and  (13),  we  may  write 

"0  {^x^y)  =  G -s{mx}  flp{u)x}@p{uJy),  (29) 

0  G~g(ujy')  Pp(u}x)  f3p(ll)y') ,  (30) 

where  G(u>)  is  given  by  (22)  for  d  G  {1,  2}.  Choosing 

X  =  Ks(ux)Ti{uy^Pp(uJx)Pp(u)y^,  (31) 

X  (wx,  t-Oy)  Ks  {u)y^jT\  (oJx)(dp {iOx)fdp (^y) ,  (32) 

5As  in  Section  2.2.2,  we  put  the  finest  scale  of  the  transform  at  m  =  0. 
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Table  5:  Impulse  responses  t\{n)  for  p  €  {0, 1, 2}. 


n 

p=0 

p=l 

p=2 

-3 

0.0078125 

-2 

0.03125 

0.046875 

-1 

0.125 

0.125 

0.1171875 

0 

0.75 

0.6875 

0.65625 

1 

0.125 

0.125 

0.1171875 

2 

0.03125 

0.046875 

3 

0.0078125 

where  K(u )  and  Ti(w)  are  digital  filter  frequency  responses,  we  may  compute  (26)  for  the 
finest  two  scales  by 

2 

^■0l(2wI,2o»j/)  x*(2w*,2 uy)  =  \(j>(ujx,ujy)\2  -  \<j>(2ux,  2ujy)\2.  (33) 

i= 1 

Inserting  the  terms  defined  by  (29),  (30),  (31),  (32),  (12),  and  (16)  with 
(f>(u)x,uy)  =  Pp(ux)PP(ujy)  into  (33)  results  in 


K{ux)G(ux)T\(uy)  +  K{uy)G{uy)Tx{ux)  +  \H{ux)\2\H{wy)\2  =  1.  (34) 


Equation  (34)  represents  a  relation  between  the  frequency  responses  of  the  digital  filters 
used  to  implement  a  multidimensional  discrete  dyadic  wavelet  transform  and  is  a 
multidimensional  analog  to  (17). 

Solving  (34)  for  Tx{uj)  by  substituting  K(uj)G(u>)  from  (17)  yields  the  closed  formula  [11] 

TiM  =  2  (i  +  I^MI2)  (35) 

In  Table  5  we  provide  the  filter  coefficients  for  Ti(oj)  from  (35)  for  p  G  {0, 1, 2}.  All  other 
filters  from  (34)  were  already  specified  in  Section  2.2.2. 

As  in  the  one-dimensional  case,  a  two-dimensional  discrete  dyadic  wavelet  transform  can  be 
implemented  as  a  fast  hierarchical  filtering  scheme.  To  derive  such  an  implementation,  we, 
similar  to  the  one-dimensional  case  from  Section  2.2.2,  use  the  definition  of  the 
two-dimensional  dyadic  wavelet  transform  (25)  and  require  s(u)x,ujy )  =  0  for  \ux\  >  n  or 
\uy\  >  7 r.  Using  Shannon’s  sampling  theorem  in  two  dimensions  [38],  (29),  (30),  and  m= 0, 
we  get 


Wr 


POO  POi 

d  *(*>  y)=  / 

J  —03  J  — C 


oo  oo 

H  YL  s(*z>  h)  sinc(4  -  ix)  Sine (ty  -  iy)- 

Zjj  —  OO  Xy  —  OO 


oo 

•  Yj  9-s(m)pp(x -tx-m)Pp(y -ty)dtxdty 

m~ — oo 
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and 


W, 


POO  p  cx 

o  s(*»  y)=  / 

«/—  00  J—  c 


OO  CX) 

X]  53  s(**>  ®»)  sinc(4  -  **)  sinc^  -  ij,)- 

lX“  —  OO  Zy  —  —  CXD 


OO 

‘ tj,)  ^  ]  Q—sijri) Ppiy  Zy  m.)  dtx  dtp. 

m=— oo 


As  in  one  dimension,  we  make  use  of  the  fact  that  the  cardinal  spline  functions  converge  to 
the  sine  function  as  their  order  r  tends  to  infinity,  and  write 


Wq  s(ujx,u)y)  —  S(^u}x,ojy)Br  (ujx)Br  (wy )/3p_|_r_j_i (wx)/3p_j_r_|_i G —s (wx) 


and 

Wq  s(u)x,  UJy)  —  S(oJX,U]y)Br  (u>x)Bt  (u>y')Pp-\.r  +  l(u!x)Pp+r+l(u)y)  G -.s(u)y)  , 

or  for  m  G  [0,  M)  and  discrete  signal  processing 

T{W^s{x, y) |  }  S(ux,  ujy)  B~l(ux)  B;l{ujy)  Bp+r+l(ux)- 

IX— -Tlx  )2/ - ZZ y 

m— 1 

•■Sp+r+i(wy)  G-S(2mux)  n  H-'iVufJH-'pUy).  (36) 

i=0 

and 

F{  WnS{x,  y)|  _  }  ^  S{ux,uy)  B;l{ux)  B~l(wy)  Bp+r+1(u>x)- 

I X — 77.  x  jZ/ — ZZy 

m— 1 

•Sp+r+iK)  C?_s(2mu;y)  ]J  (37) 

i=0 

Equations  (36)  and  (37)  describe  the  decomposition  part  of  the  filter  bank  implementation 
of  a  two-dimensional  discrete  dyadic  wavelet  transform.  The  reconstruction  part  can  be 
obtained  from  (31)-(33)  with  (12)  and  (16).  The  entire  filter  bank  implementation  of  the 
transform  is  shown  in  Figure  9.  Except  for  the  prefiltering  and  postfiltering,  we  readily 
recognize  the  implementation  proposed  in  [11]. 

Using  the  fact  that  a  wavelet  tp(x)  is  equal  to  a  first  (d=  1)  or  a  second  (d= 2)  derivative  of 
a  spline  function  fip+d(x),  (27)  and  (28)  may  be  rewritten  as 

^V)  =  -i Me  {1,2}, 

where 

=  PP+d(x)  (3P{y) 

and 

d2(x,y)  =pp{x)  Pp+d(y)- 
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(b) 


Figure  9:  Filter  bank  implementation  of  a  two-dimensional  discrete  dyadic  wavelet  transform 
(a)  decomposition  and  (b)  reconstruction  for  two  levels  of  analysis.  H*a( u>)  denotes  the 
complex  conjugate  of  H^s(u). 


Let  us  denote  Wms{x,y)  =  (W^s(x,y),Wls(x,y)),  V  =  (£,  ^),  A  =  V2  =  (^  +  ^), 
and  assume  that  fll(x,  y , )  can  be  approximated  by  d(x,  y)  for  both  i  £  {1, 2}. 

For  d=l  it  then  follows  that 


Wms(x,y)  =  2mV(s*'dm)(x,y).  (38) 

Thus  for  d=  2  we  can  write 

2  y)  =  22mA(s  *  i ?m)(x,  ?/).  (39) 

With  $(x,  y)  being  a  Gaussian,  finding  local  extrema  of  (38)  in  the  direction  of  gradient  V 
corresponds  to  the  filtering  stage  of  a  Canny  edge  detector  [33],  and  finding  zero-crossings 
of  (39)  corresponds  to  the  filtering  carried  out  with  a  Marr-Hildreth  edge  detector 
(Laplacian  of  Gaussian)  [39].  (Note  that  both  edge  detectors  involve  postprocessing).  Edge 
detection  based  on  finding  local  extrema  of  Wms(x,  y)  or  zero-crossings  of  £2=1  W^s(a;,  y) 
is  therefore  an  approximation  to  the  Canny  or  the  Marr-Hildreth  edge  detector  over  a 
range  of  dyadic  scales.  The  differences  stem  from  the  fact  that  fl(x,  y)  is  neither  a  Gaussian 
nor  is  i9*(rr,  y)  equal  to  fi(x,y). 

2.2.5  Steerable  Functions 

When  extending  the  transform  from  Section  2.2.2  to  two  dimensions,  one  of  the  first 
questions  that  come  to  mind  is  how  to  deal  with  the  fact  that  derivatives  can  be  defined  in 
any  direction  of  the  plane.6  In  case  of  a  first  derivative  of  a  Gaussian,  one  can  simply 
compute  first  derivatives  of  a  Gaussian  in  x  and  y  directions  and  then  determine  the 
derivative  in  any  direction  from  these  two  directional  derivatives  [33].  For  higher  order 
derivatives  of  a  Gaussian,  Freeman  and  Adelson  [34]  showed  that  order+1  directional 
operators  are  needed  for  synthesizing  the  operator  at  any  orientation.  In  fact,  functions 
with  the  property  of  expressing  their  arbitrary  rotations  as  linear  combinations  of  some 
functions  are  not  limited  to  derivatives  of  a  Gaussian.  Below,  we  briefly  restate  some  of  the 
results  from  [34]. 

A  two-dimensional  function  is  called  “steerable”  if  its  rotations  generate  a  finite 
dimensional  space.  Rotations  of  a  steerable  function  /(r,  0)  can  therefore  be  expressed  as 

/M-0o)  =  i>(0o)/iM ).  (40) 

i~l 

6Second  derivatives  of  central  B-splines  can  be  used,  as  we  saw  in  Section  2.2.4,  to  approximate  Laplacian 
of  a  Gaussian  for  approximately  rotation-invariant,  although  not  directional,  processing. 
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where  90  denotes  an  arbitrary  angle,  {cj(0o)}  stands  for  a  set  of  interpolating  functions, 
{/i(r,  #)}  is  a  set  of  basis  functions,  and  r  =  ^/x1  +  y1  and  6  =  arg(a;,  y)  are  polar  radius 
and  angle,  respectively. 

If  / (r,  6)  represents  a  filter  kernel,  the  result  of  filtering  with  a  rotated  filter  f(r,  9  —  0O)  can 
be  computed  simply  by  {cj(0o)}  weighted  linear  combination  of  outputs  from  basis  filters 
{/i(r>0)}-  Only  the  outputs  from  basis  filters  need  to  be  precomputed  and  then  the  output 
from  a  filter  rotated  by  any  angle  90  can  be  found  by  interpolating  between  them.  When  a 
large  number  of  rotations  of  a  template  filter  is  required,  the  above  scheme  can  lead  to 
substantial  savings  in  both  computational  cost  (time)  and  memory  requirements  (space). 
The  necessary  condition  for  a  function  f(r,9 )  to  be  steerable  is  that  f(r,9 )  is  bandlimited 
in  its  polar  angle: 

/M)=  £  «»(r)e'"*.  (41) 

n=-N 

This  can  be  verified  by  inserting  (41)  into  (40)  and  by  assuming,  for  convenience,  that 
fi(r,  9)  =  f(r,9-9i ): 

an{r)  e~jn0°  =  ]T  Ci(90)an(r )  e~3n6i,  (42) 

i= 1 

where  {#*}  is  a  set  of  user  defined  angles  and  n  G  [-N,  0].  7  Since  only  nonzero  coefficients 
an(r)  are  of  interest,  both  sides  of  (42)  can  be  divided  by  an(r).  This  yields  a  matrix 
equation  from  which  interpolating  functions  Ci(90 )  can  be  determined: 


1  ' 

’  1  1  •••  1  " 

ci(90) 

ej9° 

— 

eJ0i  2  . . .  ej0i 

C2(0o) 

JN0o 

ejN0i  ^3  NO 2  .  .  ,  gjNOi 

.  C/(0O)  . 

(43) 


For  coefficients  an  =  0  the  rows  corresponding  to  each  n  were  removed  from  the  matrix 
formulation  shown  in  (43).  For  this  system  to  have  a  solution,  the  angles  {0*}  must  be 
chosen  such  that  the  columns  of  the  matrix  are  linearly  independent. 

In  [34]  they  proved  that  the  minimum  number  of  basis  functions  fi(r,9)  needed  to  steer 
f{r,9)  according  to  (40)  is  equal  to  the  number  of  nonzero  coefficients  an(r)  in  the  Fourier 
series  expansion  (41). 

To  conclude  this  brief  description  of  steerability,  let  us  only  remark  that  functions  which 
are  not  steerable  (i.e.,  do  not  have  a  finite  number  of  terms  in  (41))  can  be  approximated 
with  steerable  functions  (a  singular  value  decomposition  was  employed  for  approximating 
such  functions  efficiently  [40]),  and  that  the  technique  of  expressing  transformed  versions  of 

7Note  that  the  constraints  are  the  same  for  n  €  [-N, -1]  and  n  €  [1,1V],  so  that  a  subset  of  all  possible 
values  for  n  6  [ —N ,  N]  can  be  taken. 
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a  function  as  linear  combinations  of  a  fixed  set  of  basis  functions  is  not  limited  to  rotations 
(extensions  to  translations  [41],  scalings  [40,  41],  and  general  transformations  [42]  have 
been  reported). 

2.2.6  Steerable  Dyadic  Wavelet  Transform 

Returning  to  the  question  from  the  beginning  of  Section  2.2.5,  the  answer  seems  obvious: 
one  needs  to  construct  a  steerable  analog  to  the  one-dimensional  transform  from  Section 
2.2.2.  Steerable  transforms  are  nothing  new — quite  a  few  [43,  44,  45,  41]  have  been  devised, 
some  of  them  [44,  45]  exactly  for  the  computation  of  directional  derivatives.  Here,  we  are 
not  interested  in  any  directional  derivatives:  we  want  to  use  derivatives  of  central  B-splines 
which,  as  the  order  of  B-splines  increases,  tend  to  derivatives  of  a  Gaussian. 

We  define  a  steerable  dyadic  wavelet  transform  of  a  function  s(x,  y)  £  L2(R 2)  at  a  scale 
2m,  m  £  Z,  as  [9] 

fOO*,  y)  =  s  *  y ),  (44) 

where  tplm(x,y)  denotes  ipm{x,y)  rotated  by  9U  ipm{x,y)  =  2~2mijj(2~mx,2~my))  ip(x,y)  is  a 
steerable  wavelet  that  can  be  steered  with  I  basis  functions,  and  9i  =  with 
ie  {1,2,...,/}. 

Analogously  to  the  one-dimensional  case  (cf.  Section  2.2.2)  we  require  the  two-dimensional 
Fourier  plane  to  be  covered  by  the  dyadic  dilations  of  'tpl{2mujx,  2muy):  there  must  exist 
A3  >  0  and  B3  <  oo  such  that 

OO  I 

<  E  B,  (45) 

m=— oo  i=l 

is  satisfied  almost  everywhere. 

If  (nonunique)  reconstructing  functions  Xm(x>  v)  are  chosen  such  that  their  Fourier 
transforms  satisfy 

OO  I 

£  =  (46) 

m=— oo  i= 1 

the  function  s(x,  y)  may  be  reconstructed  from  its  steerable  dyadic  wavelet  transform  by 

OO  / 

s(x,y)=  'HwLs*xtn(x,y),  (47) 

771— —  OO  i—1 

where  xL(x,y)  denotes  Xm(x,y)  rotated  by  0*  and  Xm(x,y)  =  2-2mx(2~mx,2 ~my). 

To  derive  an  algorithm  for  the  fast  computation  of  the  transform,  we,  similar  to  (10), 
introduce  two  smoothing  functions  such  that 

OO  / 

0K,  uy)  0(ujx,  uy)=Y,Y,  2 mu>y)  xi(2mu;1, 2muy).  (48) 

m= 0  i= 1 
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We  choose  wavelets  that  are  steerable  analogs  to  the  one-dimensional  wavelets  from  Section 

2.2.2:8 

„  / cinf  iiir. \  \ 

Ip{ur,uj0)  =  (jur  cos(ue))  y-^-J  ,  (49) 

where  ur  =  +  c u%  and  ug  =  arg (cox,u)y).  These  wavelets  can  be  steered  with  d+1  basis 

functions  (i.e.,  I  in  (40)  is  equal  to  d+1). 


Choosing 

4>(2ur)  =  Hst(u}r)  (t>{ur), 

(50) 

1pl(uT,Ug)  =  G st{u)T ,  LOg  —  6j)  (j>(oJr), 

(51) 

<p(2tvr)  =  Lst(ujr )  <p(ur), 

(52) 

and 

Xl(ur,Ug)  =  Kst(0Jr,Ulg  -  6i)  <p(0Jr), 

(53) 

and  using  (50)  through  (53)  with  (48)  computed  for  the  finest  two  scales,  we  obtain 

I 

X]  G st (kV ?  ^ <b  ) Kst (u)r  j  0i)  Hst(ur)Lst(ur)  =  1. 

i=  1 

(54) 

Setting  < j>(ur )  =  /3p(ur),  and  employing  (12)  and  (49)  with  (50)  and  (51),  we  find  that 


Hstif^r)  —  H—a{u)r) 

and 

Gst(ujr,u}g )  =  (cos(uje))dG-s(u}r), 

where  H(u)  and  G(u>)  are  specified  by  (16)  and  (22),  respectively. 

By  inserting  (55)  and  (56)  into  (54),  the  missing  two  filters  can  be  chosen  as 

Lst  (<+)  =  Ls(u>r ) 


(55) 

(56) 


(57) 


and 


Kst(ur,Ug )  =  —(cOs(wg))dKs(u)r), 


where  L(u)  and  K(u)  are  given  by  (23)  and  (24),  respectively,  and  C<i  = 

E'=,(<™(^ -»,))“ 


(58) 


Figure  10  shows  a  mammogram  and  the  corresponding  wavelet  transform  coefficients  for 
the  second  derivative  (d=2)  wavelet.  Wavelet  transform  coefficients  for  the  fourth 
derivative  (d=4)  wavelet  are  shown  in  Figure  11.  With  increasing  order  of  the  derivative, 
the  orientation  selectivity  becomes  better. 


8This  choice  can  be  viewed  as  an  extension  of  the  transform  presented  in  [46,  9,  47]. 
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Figure  11:  Wavelet  coefficients  for  mammogram  from  Figure  10(a).  d= 4,  m  €  {2,3,4}  (left 
to  right),  and  0*  =  h^7t,  i  €  (1, 2, 3, 4, 5}  (top  to  bottom). 
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2.2.7  Multiscale  Spline  Derivative-Based  Transform 

Let  us  pause  here  for  a  brief  assessment  of  the  two-dimensional  steerable  transform  derived 
so  far.  We  have  chosen  steerable  wavelets  (49)  which  are  equal  to  d-th  order  derivatives  of 
circularly  symmetric  spline  functions  in  the  direction  of  x-axis  (note  that  knots  for  these 
splines  are  circles)  and  we  have  laid  a  foundation  for  filter  bank  implementations  in  (54). 
An  obvious  shortcoming  of  this  scheme  is  the  fact  that  none  of  the  filter  kernels  obtained 
from  (55)  through  (58)  is  compactly  supported  on  the  rectangular  grid.  For 
implementations  using  digital  filters,  one  is  therefore  forced  to  approximate  these  frequency 
responses,  and  by  doing  so,  (54)  may  not  hold  anymore.  Filters  in  filter  bank 
implementations  of  steerable  pyramids  described  in  [44,  45,  41],  for  example,  were  designed 
by  using  various  techniques  for  approximating  desired  frequency  responses.  None  of  the 
reported  filter  banks  achieved  perfect  reconstruction  and  all  filter  kernels  were 
nonseparable.  Here,  we  will  take  a  different  approach.  We  will  approximate  the  wavelets  in 
(49)  in  a  way  that  will  lead  to  x-y  separable  filters  in  a  perfect  reconstruction  filter  bank 
implementation  of  the  transform  such  that  the  quality  of  approximation  will  improve  by 
increasing  the  order  of  B-splines. 

Let  us  approximate  wavelets  from  (49)  with 

r^-Pp+div)-  (59) 

Based  on  the  fact  that  B-splines  tend  to  a  Gaussian  as  their  order  increases,  it  is  easy  to 
see  that  both  wavelets  (49)  and  (59)  converge  to  the  same  functions  (i.e.,  d-th  order 
derivatives  of  the  normalized  Gaussian  in  the  direction  of  rr-axis)  as  p  — >  oo. 

In  order  to  steer  wavelets  ip(x,  y)  given  by  (59)  (note  that  steering  will  be  only 
approximate,  since  these  wavelets  are  not  steerable),  we  need  to  find  basis  functions  that 
will  approximately  steer  ip(x,y).  Computing  rotations,  as  we  did  in  (44),  is  not  practical 
here,  because  arbitrary  rotations  of  (59)  cannot  be  expressed  exactly  in  terms  of  central 
B-spline  functions  from  Section  2.1.  Instead,  we  take  advantage  of  the  property  of 
circularly  symmetric  functions  that  rotations  of  their  directional  derivatives  are  equal  to 
directional  derivatives  in  rotated  directions: 

f  d^Qcjx,  y) }  _  d^Qcjx,  y) 

00  \  dnd  J  dnd$o  ’ 

where  TZg0  stands  for  rotation  by  9q,  deCQ^  =  n'Vgc(x, y),  gc(x,y)  is  a  circularly  symmetric 
function,  and  fig0  denotes  vector  n  =  (cos  9,  sin  6)  rotated  by  9q. 

Let  us  choose 

Q(x,y)  =  Pp+d(x)flp+d{y), 
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which  is  approximately  circularly  symmetric  function  for  higher  order  splines.  A  rotation  of 
ip(x,  y)  =  from  (59)  by  9q  can  therefore  be  approximated  by 


_  =  y,  (  d  \  d_j  ,  dd  i/3p+d(x)  d^p+d{y) 

dnd  £^\i)x  y  dx d~i  dy * 


where  n  =  (cos0o>sin0o)  =  {nx,ny). 

Note  that  in  case  of  Gaussian,  which  is  both  x-y  separable  and  circularly  symmetric,  (60) 
becomes  exact  (e.g.,  for  g(x,  y)  =  e~(x2+y2\  90  =  -9,  and  d  =  {2, 4},  we  obtain,  up  to  a 
scaling  factor,  x-y  separable  basis  set  for  the  second  and  fourth  derivative  of  Gaussian  from 
Tables  III  and  VII  of  [34]). 

Having  found  a  set  of  basis  functions  (60)  that  approximately  steer  wavelets  (59),  we  want 
to  construct  a  transform  such  that  Equations  (44)  through  (48)  will  be  valid  (superscript  i 
must  be  viewed  now  as  an  index,  rather  than  rotation  by  0*).  In  frequency  domain,  we  can 
express  basis  functions  from  (60)  as 

i’t+1(vx,Wy)  =  GdTst(ojx)Gl_s(ujy)Pp+i(ujx)Pp+d-i(uJy),  i  =  0, 1, . . . ,  d,  (61) 


where  Gd(u)  is  given  by  (22)  and  G°(u)  =  1. 

Choosing  appropriate  xj(cjx,  vy)  to  obtain  a  relation  needed  for  the  filter  bank 
implementation  of  the  transform  is  more  complicated  than  in  one  dimension.  Since  we 
could  not  find  a  general  solution  for  arbitrary  d ,  we  solve  each  case  separately.  Below,  we 
present  solutions  for  the  first  four  orders.  When  d  <  2,  we  impose 

<p(ujx,u)y)  =  (j)(oox,u)y)  =  Pp(ux)Pp(LOy),  a  constraint  analogous  to  the  one  from  Section  2.2.2. 
For  d  =  1,  we  write  similar  to  [11] 

X  =  dCs{uj^Tii^y^^p{u)xsjpp—i{u}y),  (62) 

X  =  Ti(u)x)Ks(LOy)Pp—i(u)x){3p(u!y),  (63) 

where  Kd(u>)  and  Ti(cj)  are  given  by  (24)  and  (35),  respectively. 

Computing  (48)  for  the  finest  two  scales  and  inserting  (12),  (61),  (62),  and  (63)  yields  a 
relation  between  frequency  responses 


GHoo^KUuj^T ;  U.)  +  T,  fuUG1! 


For  d  =  2,  we  choose 


X  (wx,  LOy)  -  Ks(u)x)T2(u>y)Pp(u>x)j3p—2(k>y), 

X2(uJX,UJy)  =  Kl(Ux)Kl(uJy)Pp-l(U}X)Pp-l(Uy), 

X  {^X  >  Uy)  =  T2  (u)X  )  Ks  ( U)y  )  ftp- 2  i^X  )  Pp  {k)y  )  , 
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(64) 

(65) 

(66) 


where 


(67) 


T2(u)  =  \H(u)\\ 

Using  (12),  (61),  and  (64)  through  (66)  with  (48)  results  in 

G2  (ux)K2  (u}x)T2  (u>y)  +  G'MK'MG'MK'M  +  T2(ux)G2(cvy)K2(uy)+ 

+\H(ujx)H(uy)\2  =  1. 

For  orders  d  >  2,  we  require  (j)(uix,u)y)  =  /3p(cux)  j3p(uy)  and  (p(tox,uy)  =  (p(ujx)cp(ujy),  where 


<p(u>)  is  specified  by  (14)  and  (23). 

For  d  =  3,  we  choose  reconstructing  functions 

Xl(wx,wy)  =  KsMvMv- 3(w„),  (68) 

X2(<jjx,uy)  =  -K](u}x)Kl(uy)V3(ivx)V3(uy)<p-i(uJx)0-2M,  (69) 

X3{Vx,Vy)  =  -Kl(ux)K2(uy)V3{ux)V3(u>y)<P-2{oJx)<P-l(Uy),  (70) 

x4(ux,uy)  =  K3s(u}y)0-3{ux)ip(ojy),  (71) 

where 

UM  =  ^(1  -  lifMP),  (72) 


and  (p-i(uj )  G  Ll(R)  denotes  a  function  such  that  <p(u>)  =  /3i_i(u>)<£_i(u>),  i  G  N. 
Employing  (61),  (12),  (14),  and  (68)  through  (71)  with  (48)  yields  a  relation 

G3K)tf3K)  -  G2(ux)K2(ux)V3(ux)G1(u}y)K1(u)y)V3(ujy)— 
—G1(u>x)Kl(u}x)V3(ujx)G2(u)y)K2(ujy)V3(u>y)  +  G3(coy)K3(ujy)+ 


~\~H  (u)x}Ij(u)x')H  (u)y}L{u)y'j  —  1, 

where  L(u)  is  specified  by  (23). 

For  d  =  4,  our  choices  are 

xH^x^y)  =  K4S  (ujx)T2(ujy)(f(u}x)<p^  (ojy ) ,  (73) 

X2(tvx,uy)  =  K*(ux)Kl(u)y)<p-i(<ox)$-3(u>y),  (74) 

X3(u)x,uy)  =  —K2  (lox)K2  (uy)  V4  {ujx)V4  {ujy)(p-2  (ux)<P-2  (uy) ,  (75) 

X4{u>x,Uy)  =  Kl(ux)K*(wy)<p-a(wx)0-i(u)y)t  (76) 

fk.Wj,)  =  T2{ux)K*(ujy)(p-i(vx)0(ujy),  (77) 

where 

VA{u>x)  =  l-\H(u>)\2.  (78) 
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Using  the  above  functions  (73)  through  (77),  (61),  (12),  and  (14)  in  (48)  computed  for  the 
finest  two  scales  gives 


G4  {lox)Ka  (wx)T2  (uy)  +  G3  (<jjx  )K3(u>x)G1(uy)K1(uy)— 
—G2  {uix)K2  (ux)V±  (ujx  )  G2  (ojy )  K2  (uy )  U4  (u)y )  +  G1  (u>x  )K1(ux)G3(uy)K3  (uy ) + 

+T2{ux)G\uy)K\uy)  +  H(ux)L(ux)H(uy)L(uy)  =  1. 


Here,  we  have  even  more  freedom  for  choosing  the  reconstructing  functions  than  in  one 
dimension.  The  above  functions  for  d  =  {2,3,4}  were  found  by  trying  to  imitate  the 
one-dimensional  transform  from  Section  2.2.2  as  much  as  possible.  All  decomposition  filters 
Gd(ui)  were  first  paired  with  corresponding  reconstruction  filters  Kd(u),  and  then  all  other 
potential  digital  filters  were  specified  as  polynomials  in  i7_s(w).  We  inserted  thus  specified 
filters  into  a  relation  between  their  frequency  responses  and  solved  for  the  unknown 
polynomial  coefficients.  Since  we  allowed  more  filters  with  higher-degree  polynomials  than 
expected  in  the  solution,  the  resulting  system  of  linear  equations  was  underdetermined. 
This  allowed  enough  freedom  for  removal  of  undesired  digital  filters  and  for  balance 
between  degrees  of  polynomials. 

The  described  procedure  for  determination  of  reconstructing  functions  and  filters  involves 
quite  a  lot  of  heuristics  to  obtain  the  appropriate  solution  from  the  underdetermined  linear 
system.  Unfortunately,  we  are  not  aware  of  any  systematic  way  (aside  from  numerical 
optimization,  which  may  be  pretty  cumbersome)  to  obtain  solutions  comparable  to  the 
ones  above. 

Next,  we  will  derive  a  filter  bank  implementation  of  the  transform.  As  in  Section  2.2.4,  we 
assume  a  bandlimited  input  signal:  s(ux,u)y)  =  0  for  \ux\  >  n  or  \u)y\  >  7r.  Using  Shannon’s 
sampling  theorem  in  two  dimensions  [38]  with  (44)  and  basis  functions  from  (61),  we  can 
write 


W, 


r  oo  roo  w  w 

o®(®,  V)~  /  X)  X  s(**’  h)  sinc(^  -  **)  sinc(£j,  -  iy)- 

J — oo  J  — OO  *  ■  


%x —  OO  Zy —  OO 


oo  oo 

^  /  9— s ' TP'x)  O—siV^y^fip+d—iiy  £y  ^h/)  dtx  dty , 

mx=—oo  my=— oo 


where  i  =  0, 1, . . . ,  d  as  in  (61). 

Again,  we  approximate  sine  functions  with  r-order  cardinal  splines,  then  use  (4),  and  get 


DFT{  W^s(x,  y) \  }  ~  S(ux,u>y)  Br  l{ux)  Br  V^)  Bp+r+i+l(ux)- 

I  JC — Tlx  jj/ — »£ y 

m— 1 

i+iK)  Gl7(2mcua:)  GU^y)  II  Hp_+si(2ncvx)Hp_+sd-i(2nUy).  (79) 
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Using  (79)  with  an  approximation  Bp+r+i+l( oS)  ~  Bp+r(co)Bi(u),  we  can  obtain  a  filter 
bank  implementation  of  the  transform  decomposition.  The  reconstruction  part  follows  from 
(48),  (61),  and  reconstructing  functions  for  distinct  values  of  d.  Figure  12  shows  filter  bank 
implementations  of  a  multiscale  spline  derivative-based  transform  for  d  =  {1, 2, 3, 4}.  For 
d  =  1,  we  recognize  (except  for  the  prefiltering  and  postfiltering)  the  filter  bank 
implementation  of  a  two-dimensional  discrete  dyadic  wavelet  transform  from  [11].  For 
d  =  2,  however,  our  transform  differs  from  the  filter  bank  presented  in  [8]  (i.e.,  the 
corresponding  transform  described  in  Section  2.2.4):  second  derivative  is  computed  only  in 
the  directions  of  x  and  y- axis  in  [12,  8],  which  is  not  enough  for  steering.  Although  not 
particularly  appropriate  for  orientation  analysis,  such  a  scheme  may  still,  as  we  have  seen 
in  Section  2.2.4,  efficiently  approximate  Laplacian  of  Gaussian  across  dyadic  scales. 

A  transform  similar  to  the  one  described  in  this  section,  was  presented  in  [44,  45,  41]. 

Their  filter  bank  implementation  is  less  redundant  (downsampling  is  used  on  the  lowpass 
branch,  while  simultaneously  keeping  aliasing  negligible  by  using  a  filter  with  insignificant 
amount  of  energy  for  |uv|  >  f)  and  uses  reconstruction  filters  with  same  magnitude 
frequency  responses  as  the  decomposition  ones — a  possible  advantage  for  certain 
applications.  They  have,  on  the  other  hand,  little  control  over  the  function  from  which 
derivatives  are  computed  (to  obtain  a  d-th  derivative,  they  multiply  a  circularly  symmetric 
filter  by  (—j  cos  9)d  with  all  filters  being  obtained  by  recursive  minimization  of  a  weighted 
combination  of  constraints),  filter  bank  does  not  have  perfect  reconstruction,  and  none  of 
the  filters  is  x-y  separable. 

2.2.8  Finite  Impulse  Response  Filters 

Since  all  two-dimensional  filters  used  in  the  filter  bank  implementations  of  the  transforms 
from  previous  sections  are  x-y  separable,  we  will  begin  this  section  with  a  detailed 
description  of  FIR  filter  implementations  for  the  one-dimensional  discrete  dyadic  wavelet 
transform  implementation  from  Figure  6.  The  extension  to  two  dimensions  will  then  be 
straightforward . 

Let  us  refer  to  filters  applied  at  scale  2TO  as  filters  at  level  ra+1,  and  let  filters  at  level  1 
(Equations  (16),  (22)  through  (24),  (35),  (67),  (72),  and  (78))  be  called  “original  filters,”  to 
distinguish  them  from  their  upsampled  versions.  As  an  input  to  the  filter  bank  from  Figure 
6,  we  consider  a  real  signal  s(n)  e  l2(Z),  n  e  [0,  iV  -  1],  Depending  on  the  length  of  each 
filter  impulse  response,  filtering  an  input  signal  may  be  computed  either  by  multiplying  the 
discrete  Fourier  transforms  of  the  two  sequences  or  by  circularly  convolving  s(n)  with  a 
filter’s  impulse  response.9  Of  course,  such  a  periodically  extended  signal  may  change 

9As  is  customary  in  image  processing,  we  use  circular,  rather  than  linear,  convolution. 
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b;K)  B>y)  - 

-  Bp+^og  Bp+r+1(coy) 

(a) 

Bp+r+l(®x)  Bp+r+i(©y) 

it  Br(cox)  Br(wy) 

(b) 


(c)  (d) 


(e)  (f) 


Figure  12:  Filter  bank  implementation  of  a  multiscale  spline  derivative-based  transform  for 
[0,M  —  1]:  (a)  Prefiltering,  (b)  postfiltering,  (c)  decomposition  and  (d)  reconstruction 
modules  for  d  =  1,  and  (e)  decomposition  and  (f)  reconstruction  modules  for  d  =  2. 
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Figure  12:  Continued:  (g)  Decomposition  and  (h)  reconstruction  modules  for  d  —  3,  and 
(i)  decomposition  and  (j)  reconstruction  modules  for  d  =  4.  Decomposition  modules  are 
recursively  applied  at  the  locations  of  the  filled  circles. 
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abruptly  at  the  boundaries  causing  artifacts.  A  common  remedy  for  such  a  problem  is 
realized  by  constructing  a  mirror  extended  signal  [12] 


s(— n  —  1)  if  n  €  [— N,—  1] 
s(n)  if  n  G  [0,  N  —  1], 


(80) 


where  we  chose  the  signal  sme(n)  to  be  supported  in  [-N,  N  -  1].  It  will  become  evident 
shortly,  that  mirror  extension  is  particularly  elegant  in  conjunction  with 
symmetric/antisymmetric  filters. 

Let  us  first  classify  symmetric/antisymmetric  real  even-length  signals  into  four  types  [23]: 


Type  I  f(n)  =  /(-n), 


Type  II  f(n)  =  f(—n  —  1), 


Type  III  f(n)  =  -f(-n), 

Type  IV  f(n)  =  — /(— n  —  1), 

where  n  G  [— N,  N  —  1].  Note  that  for  Type  I  signals  the  values  at  /( 0)  and  f(—N )  are 
unique,  and  that  for  Type  III  signals  the  values  at  /( 0)  and  /(—AT)  are  equal  to  zero. 

Using  properties  of  the  Fourier  transform  it  is  easy  to  show  that  the  convolution  of 
symmetric/antisymmetric  real  signals  results  in  a  symmetric/antisymmetric  real  signal.  If 
a  symmetric/antisymmetric  real  signal  has  an  even  length,  then  there  always  exists  an 
integer  shift  such  that  the  shifted  signal  belongs  to  one  of  the  above  types. 

Now,  we  are  ready  to  examine  the  filter  bank  implementation  of  a  one-dimensional  discrete 
dyadic  wavelet  transform  from  Figure  6  with  filters  given  by  (16)  and  (22)  through  (24) 
driven  by  a  mirrored  signal  sme(n)  at  the  input.  Let  the  number  of  levels  M  be  restricted 

by 

M  <  l+log^-*  (81) 

where  Lmax  is  the  length  of  the  longest  original  FIR  filter  impulse  response. 

Each  FIR  filter  block  in  the  filter  bank  consists  of  a  filter  and  a  circular  shift  operator. 
Equation  (81)  guarantees  that  the  length  of  the  filter  impulse  response  does  not  exceed  the 
length  of  the  signal  at  any  block. 

Since  our  input  signal  sme(n )  is  of  Type  II  and  noninteger  shifts  at  level  1  are  rounded  to 
the  nearest  integer,  it  follows  that  a  processed  signal  at  any  point  in  the  filter  bank  belongs 
to  one  of  the  types  defined  above.  This  means  that  filtering  a  signal  of  length  2 N  can  be 
reduced  to  filtering  a  signal  of  approximately  one  half  of  its  length.  (For  Types  I  and  III, 

N  +  1  samples  are  needed.  However,  for  Type  III  one  needs  to  store  only  N  —  1  values 
because  zero  values  are  always  present  at  the  zeroth  and  (— iV)-th  sample  position). 
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Implementation  is  particularly  simple  for  FIR  filters  designed  with  d  even  and  p  odd. 
Filters  are  of  Type  I  in  this  case,  so  the  signal  at  any  point  of  the  filter  bank  will  be  of 
Type  II.  An  FIR  filter  block  from  the  filter  bank  shown  in  Figure  6  can  therefore  be 
implemented  by 


L—  1 
2 


Ts,mw(n) 

=  /(0)«//(n)  +  f(i)[un{n 

i—1 

-  2mi)  +un(n  +  2mi)],  n€[0,iV-l], 

(82) 

where 

u(—n  — 

1)  if  n  €  [— y,  — 1] 

uu(n)  =  <  u(ri) 

if  n  e  [0,  N  -  1] 

(83) 

k  u{2N- 

n  —  1)  ifne[N,*£], 

u(n)  is  an  input  signal  to  a  block,  f(n)  is  an  impulse  response  of  some  original  filter,  L  is 
the  length  of  the  filter,  and  N  is  the  length  of  an  input  signal  s(n)  to  the  filter  bank. 
Implementation  of  filters  bp(n )  used  for  prefiltering  and  postfiltering  represents  a  special 
case  of  (82)  with  m=0. 

A  filter  bank  with  the  above  implementation  of  blocks  and  signal  s(n)  at  the  input  yields 
equivalent  results  as  circular  convolution  for  sme(n )  as  defined  by  (80).  In  addition  to 
requiring  one  half  the  amount  of  memory,  the  computational  savings  over  a  circular 
convolution  implementation  of  blocks  are,  depending  on  the  original  filter  length,  three  to 
four  times  fewer  multiplications  and  one  half  as  many  additions. 

A  similar  approach  can  be  used  for  other  filters.  However,  things  get  slightly  more 
complicated  in  this  case,  because  the  filters  are  not  of  the  same  type  and  the  signal 
components  within  the  filter  bank  are  of  distinct  types.  As  a  consequence,  an 
implementation  of  blocks  that  use  distinct  original  filters  may  not  be  the  same,  and  the 
implementation  of  blocks  at  level  1  may  differ  from  the  implementation  of  blocks  at  other 
levels  of  analysis. 

The  decomposition  blocks  at  level  1  can  be  implemented  by 

f-i 

G-Sflu(n)  =  ^2  0(*)[«//(n  -  *  -  1)  -  un{n  +  <)],  n  €  [1,  N  —  1], 

i= 0 

for  d  odd,  (82)  for  d  even, 

£_i 
2  1 

H-St0u(n)  =  h(i)[un(n  -  i  -  1)  +  un(n  +  *)],  n  e  [0,  N ], 

i= 0 

for  p  even,  and  (82)  for  p  odd,  where  un{l)  is  defined  by  (83),  g(n )  and  h(n)  are  impulse 
responses  of  the  filters  computed  from  (22)  and  (16),  respectively,  and  L  is  the  length  of 
the  corresponding  impulse  response. 


The  output  from  a  block  G-S(u )  at  level  1  is  of  Type  III  for  d  odd  and  of  Type  II  for  d 
even,  while  the  output  from  H^s(u)  at  the  same  level  is  of  Type  I  for  p  even  and  of  Type  II 
for  p  odd. 

The  decomposition  blocks  at  subsequent  levels  m  G  [1,M— 1]  can  be  implemented  by 
L-\ 

G-S,mu(n)  =  g{i)[uj{n  -  2 m(i  +  s ))  -  u7(n  +  2 m(i  +  s))],  n  G  [1,N  -  1], 

i— 0 

for  d  odd  and  p  even, 

G-S,mu(n)  =  ^2  g{i)[un{n  -  2 m(i  +  s ))  -  un(n  +  2 m(i  +  s))],  n  e  [0,  N  -  1], 

i=0 

for  d  and  p  odd, 

L  —  1 

F_s>mu(n)  =  /(0)u/(n)  +  ^  /(i)[u7(n  -  2 mi)  +  it7(n  +  2TO*)],  n  G  [0,  iV],  (84) 

i=  1 

with  /(n)  =  p(n)  for  d  and  p  even, 


k_i 
2  1 


H-Stmu(n)  =  ~  2m(i  +  s))  +  u7(n  +  2m(i  +  s))],  n  G  [0,  IV],  (85) 


i=0 


for  p  even,  and  (82)  for  p  odd,  where 


ui(n )  =  < 


u(—n) 

u(n) 

u{2N  -  n ) 


if  nG  [-f,-l] 
if  n  G  [0,  AT] 
ifnG  [iV  +  1,^]. 


(86) 


The  outputs  from  blocks  G-S(2muj)  are  of  Type  III  for  d  odd  and  p  even,  of  Type  IV  for  d 
and  p  odd,  and  of  Type  I  for  d  and  p  even,  whereas  the  outputs  from  H_s(2mu)  are  of 
Type  I  for  p  even  and  of  Type  II  for  p  odd. 

Next,  the  reconstruction  blocks  at  level  1  can  be  implemented  by 


L 


2 

Ksfiu(n)  =  Y2  H’i)[uin(n  -  *  +  1)  ~  um(n  +  *)], 

Z— 1 

for  d  odd,  (82)  for  d  even, 


nG  [0, iV  —  1], 


Ls,ou(n)  =  'Y  l(i)[ui(n  -i  +  1)  +  u7(n  +  *)],  n  G  [0,  AT  -  1], 
1=1 


for  p  even,  and  (82)  for  p  odd,  where 

I  -u(-n) 


um(n )  =  { 


0 

u{n ) 
0 


ifnG 
if  n  =  0 

if  n  G  [1,  AT  —  1] 
if  n  =  V 

ZN\ 


-u(2N  —  n)  ifnG  [N +  !,*£], 


(87) 
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u/(n)  is  as  defined  by  (86)  and  k(n)  is  an  impulse  response  of  the  filter  from  (24).  Note 
that  both  outputs  from  blocks  Ks(u>)  and  Ls(u)  are  of  Type  II. 

The  reconstruction  blocks  at  subsequent  levels  can  be  implemented  by 

k- 1 
2  1 

Ks,mu(n )  =  k(i  +  l)[uIH[n  -  2 m(i  +  s ))  -  uIH(n  +  2 m(i  +  s))],  n  6  [0 ,N], 

i= 0 


for  d  odd  and  p  even,  (84)  with  f(n)  —  k(n)  for  d  and  p  even, 


k_x 
2  1 


Ks,mu(n)  =  +  1  )[uiv(n  -  2 m(i  +  s))  -  uIV(n  +  2 m(i  +  s))],  n  e  [0,  N  -  1], 

z— 0 


for  d  and  p  odd, 


Ts)m‘u(? 7.)  —  H~Simu(n), 


for  p  even,  and  (82)  for  p  odd,  where  um(l)  is  given  by  (87), 


uiv(n)  =  { 


—u(—n  —  1) 
u(n) 

—u(2N  —  n  —  1) 


if  n  6  [— y,  —1] 
if  n  €  [0,  N  —  1] 
ifn€[JV,f], 


and  H^s>mu(n)  is  specified  by  (85).  We  observe  that  the  outputs  from  blocks  Ks(2mu>)  and 
Ls(2mu),  m  e  [1,  M  —  1],  are  of  Type  I  for  p  even,  and  of  Type  II  for  p  odd. 

When  we  compare  the  above  implementation  of  blocks  to  circular  convolution  driven  by  a 
mirrored  signal  sme{n)  at  the  input,  we  observe  that  approximately  twofold  less  memory 
space,  three  to  four  times  fewer  multiplications  and  one  half  as  many  additions  are 
required.  (For  Type  I  signals  an  additional  sample  has  to  be  saved  because  two  values  are 
without  a  pair). 

Until  now,  we  have  talked  only  about  the  one-dimensional  case,  whose  filter  bank 
implementation  is  depicted  in  Figure  12.  Two-dimensional  transform  filter  bank 
implementations  (Figures  9  and  12)  are  not  only  comprised  of  x-y  separable  filters  solely, 
but  also  use  all  the  filters  from  Section  2.2.2.  Everything  presented  in  this  section  so  far  is 
therefore  directly  applicable  to  the  two-dimensional  case.  Filters  which  have  not  been 
treated  yet  (i.e.,  ti(n),  t2(n),  V3 (n),  v4(n),  and  filters  bp(n)  from  the  decomposition  modules 
of  Figure  12)  can  all  be  realized  by  (82)  for  p  odd  or  m  =  0  and  (84)  otherwise  ( f(n ) 
denotes  an  impulse  response  of  any  of  the  above  mentioned  zero-phase  filters  in  this  case).10 
The  implementation  presented  in  this  section  performs  all  operations  in  the  spatial  domain, 
however,  one  could  also  implement  the  structures  shown  in  Figures  6,  9,  and  12  with  an 
input  signal  sme{n)  (Equation  (80))  in  the  frequency  domain.  For  short  filter  impulse 

10In  case  of  filters  vz  (n)  and  v4(n),  a  slightly  more  efficient  implementation  can  be  obtained  by  precom¬ 
puting  new  filters  k  *  vz (n)  and  k  *  v4(n),  and  then  implementing  them  by  Ks<mu(n),  rn  6  [0,  M  —  1]. 
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responses,  such  as  those  given  in  Tables  2,  3  and  4,  the  spatial  implementation  described  in 
this  section  is  certainly  more  efficient.  For  long  filter  impulse  responses,  however,  filtering 
is  faster  if  implemented  in  the  frequency  domain.  Additional  details  on  alternative 
implementation  strategies  can  be  found  in  [48]. 


2.2.9  Infinite  Impulse  Response  Filters 

Implementation  of  HR  filters  b~ 1(n)  which  were  introduced  in  Section  2.1.2  is  a  bit  more 
involved  than  the  one  encountered  in  the  previous  section.  Fortunately,  the  number  of 
different  cases  is  much  smaller  here:  possible  input  to  b~l  (n)  in  filter  banks  from  Figures  6, 
9,  and  12  is  either  of  Type  II  or  of  Type  I.11  We  will  use  ideas  and  a  few  results  from  [20]. 
Let  us  first  take  a  closer  look  at  the  system  function  B~1(z).  This  function  can  be  written 
as  a  cascade  of  terms 


(1  —  az~1)(  1  —  az)J 


which  can  be  expressed  in  a  parallel  form  as 


1  -  a2 


“0- 


where  a  and  ^  are  poles  of  the  causal  and  the  anticausal  filter,  respectively. 

The  impulse  response  of  this  term  can  be  written  as 

«(»)  =  r^aW- 

1  —  ar 

We  choose  to  implement  E(z)  in  a  cascade  form  and  therefore  extract  the  difference 
equations  from  (88): 

c+(n )  =  u(n)  +  ac+(n  —  1)  n  =  1, 2, . . . ,  N— 1, 


and 

c(n)  =  a  (c(n  +  1)  —  c+(n))  n  =  N— 2,  N— 3, . . . ,  0,  (91) 

where  u(n )  denotes  the  input  to  the  block,  c+(n)  is  the  output  from  the  causal  part,  and 
c(n)  stands  for  the  output  from  the  block. 

To  solve  (90)  and  (91)  we  need  boundary  conditions  c+(0)  and  c(N— 1).  Let  us  begin  with 
filters  b~ 1(n)  in  filter  bank  implementations  from  Figures  6,  9,  12(a),  12(b),  and  Figures 
12(h)  and  12(j)  with  m=0.  We  derive 

o  N- 1  i+1  ,  n2N-i  io 

c+(0)  =  53  oTluiip{i)  =  u(0)  +  53  . ,  _  2AT  m(»)  -  u(°)  +  53ai+1M(i),  (92) 

i~ — oo  i=0  ^  z=0 

11  Note  that  symmetry  types  in  this  section  slightly  differ  from  those  defined  in  Section  2.2.8:  here,  mirror 
extended  signals  are  periodically  repeated,  so  that  they  stretch  from  -oo  to  oo. 
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and,  using  parallel  form  (89) 


—a 


y  1  ryN—i  1  -yTV+l+i 

^-1)  =  Tr^(c+(W-1)  +  E  — + “  „  »(<)) 


—a: 

1  -  a2 


i=o  *  -  «2Ar 
(c+(N- 1)+  2  ^"MO). 

iz=N— 1— i0 


(93) 


where 


Uiip(n) 


■l 


uu(n  mod  (2 N))  if  n  >  0 

un{-(n  +  1)  mod  (2 N))  if  n  <  0, 


/  x  _  /  u(n)  if  n  G  [0,  N  —  1] 

un(n)  -  |  u^N  _n_y  jf  n  6  [at,  2AT  -  1], 

N  is  the  length  of  an  input  signal  to  the  filter  bank,  and  z0  <  N—l  is  selected  such  that  aio 
falls  below  a  predefined  precision  threshold. 

For  HR  filters  from  Figures  12(h)  and  12(j)  with  m  e  [1,  M  —  1]  and  p  odd,  we  get 

c+(0)  =  Y  [ar^uj/pf*)]  = 


oo  TV— 1 


(  r  I  2AT(n+l)  —  i  1 

=  u(o)  +  2^  Y  j  I a  2m  I  +  a  2m  r  u(i) 

n=o  i=o  ^  11  L  luJ 


(94) 


and 


where 


OO  TV— 1 


c(N-l)  =  ^JL.(c+(N-1)  +  E  E  (95) 


n— 0  Z— 0 


= 


ax  if  x  6  Z 
0  otherwise. 


If  A?'  is  a  power  of  two  and  2m  1  <  N  ((81)  guarantees  that  the  latter  condition  is  always 
true)  (94)  becomes 


c+(o)  =  «(o)  +  E 

z^O 


r  i±ii 

,  r  2N-ii 

I  Ot  2m 

+  la  2 m 

L  J 

u  L  J  u 

„  2jV 

1  —  a^m 


u(i), 


while  (95)  can  be  written  as 


c(N-l) 


„  ,  r  ALzii  ,  r  jv+i+n 
— n  N~l  a  2m  +  \a 

-2{c+{n-  1)  +  y  1 — ij^-L - h 


1  —  of 


1=0 


2JV 

1  —  Gl  2m 


«(*)• 


Finally,  the  boundary  conditions  for  filters  bp  1  (n)  from  Figures  12(h)  and  12(j)  with 
me  [1,  M  —  1]  and  p  even  are 


U 

c+(0)  =  Y  oT^uIp{i)^  = 
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(96) 


oo  , 

=  E{[“WF].“(o)  + 

n— 0  ^ 


JV(2n+l) 

a  2m 


u{N)+ 


2JV(n+l)— t 

a  2m 


and 


where 


c(jv)  =  r^(2c+(iV)_?i(iV))’ 


ui(n) 


uip{n)  =  «/(|n|  mod  (2N)), 

u(n )  if  n  e  [0,  iV] 

u(2N  —  n)  if  n  G  [iV  +  l,2iV  -  1], 


and  c~(N)  =  c+(N )  with  c~(n)  denoting  response  of  the  anticausal  filter  from  (89)  was 
used. 

Again,  if  N  is  a  power  of  two  and  2m_1  <  N,  (96)  can  be  simplified 


„+ 


(0)  = 


u(0 )  +  ch2™u(N)  +  Eili1  { [a2?n]u  +  }  u(i) 


2 JV 

1  —  a2m 


Series  in  expressions  for  c+(0),  c(N  —  1),  and  c(N)  for  filters  from  Figures  12(h)  and  12(j) 
with  m  €  [1,  M  —  1]  can  be,  similar  to  (92)  and  (93),  truncated  according  to  the  desired 
precision. 

For  orders  p  greater  than  three,  we  implement  B~l(z)  as  a  cascade  of  terms  E(z)  with 
different  a’s.  Note  that  the  output  from  block  E(z)  is  always  of  the  same  type  as  the  input 
to  it. 


2.3  Image  Fusion 

Image  fusion  combines  particular  aspects  of  information  from  the  same  imaging  modality  or 
from  distinct  imaging  modalities  and  can  be  used  to  improve  the  reliability  of  a  particular 
computational  vision  task  or  to  provide  a  human  observer  with  a  deeper  insight  about  the 
nature  of  observed  data.  Whether  it  is  combining  different  sensors  or  extending  the 
dynamic  range  of  a  single  sensor,  the  goal  is  to  achieve  more  accurate  inferences  that  can  be 
achieved  by  a  single  sensor  or  a  single  sensor  setting.  By  fusing  together  processed  sections 
of  images,  a  combined  image  which  is  superior  to  the  sum  of  its  parts  can  be  constructed. 
We  are  developing  methods  by  which  transformed  mammograms  can  be  broken  apart  and 
fused  together  in  a  manner  which  will  increase  the  overall  mammogram  interpretability. 
The  simplest  method  of  fusing  images  is  accomplished  by  computing  their  average.  Such  a 
technique  does  combine  features  from  input  images  in  the  fused  image,  however,  the 
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contrast  of  the  original  features  can  be  significantly  reduced.  Among  more  sophisticated 
methods,  multiscale  and  multiresolution  analyses  have  become  particularly  popular. 
Different  pyramids  [49,  50]  and  wavelet-based  techniques  [51,  46,  52,  53]  have  been  applied 
to  this  problem. 

In  this  section,  we  compare  an  image  fusion  method  based  upon  the  steerable  dyadic 
wavelet  transform  with  recently  published  fusion  methods  based  upon  the  gradient 
pyramid,  the  orthogonal  wavelet  transform,  and  the  biorthogonal  wavelet  transform. 
Sections  2.3.1  and  2.3.2  introduce  the  gradient  pyramid  and  the  discrete  wavelet  transform, 
respectively.  For  a  more  thorough  treatment  of  these  transforms,  please  refer  to  references 
cited  therein.  In  Section  2.3.3,  we  examine  the  performance  of  all  four  transforms  on 
examples  with  relation  to  mammography. 


2.3.1  Gradient  Pyramid 


Gaussian  pyramid  [30]  was  used  for  construction  of  a  gradient  pyramid  used  in  [49].  Let 
the  generating  filter  kernel  for  the  Gaussian  pyramid  be 


w{nx,  ny)  =  wb*  wb(nx,ny)  =  — 

ZD  D 


1 

4 

6 
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1 

4 
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1 

where  wb(nx,ny)  denotes  the  3  by  3  binomial  filter: 


wb(nx,ny)  —  ^ 


1  2  1 
2  4  2 
1  2  1 


Level  m  €  N  of  the  Gaussian  pyramid  for  an  input  image  matrix  s(nx,ny )  is  then 

QmS(jlxi  Wy)  =  (nJ  *  ^j/))j.2> 


with 


Gos(nx,ny )  =  s(nx,ny). 


Gradient  pyramid  is  obtained  from  the  Gaussian  pyramid  as 


/)  —  dj  *  Hy)  “t"  Wb  *  Gm^iP'xt 


where 
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d2~  V2 

d3  = 
di  =  72 


0  -1 
1  0  1  ’ 


-1 

1 


,  and 


-1  0 
0  1 


An  image  is  reconstructed  from  the  gradient  pyramid  by  converting  the  pyramid  to  the 
Laplacian  and  then  to  the  Gaussian  pyramid.  The  Laplacian  pyramid  is  approximated  as 


where 


C,ms{jlx)  Tly)  —  fly)  T  W  ♦  Kyj nS(jlxi  72.j /)) 


JCms(nx,  Uy)  —  di  *  Vms(nx,  ny). 

8  i= i 


(97) 


An  approximation  to  the  Gaussian  pyramid  is  obtained  by 

Uy)  =  dlms{nx,  Tly )  +  Aw  *  Tly))^2- 


Note  that,  since  the  filters  di  have  the  center  of  symmetry  between  samples,  they  need  to 
be  shifted  for  reconstruction,  and  that,  because  of  the  approximation  (97),  the  gradient 
pyramid  does  not  have  the  perfect  reconstruction  property. 


2.3.2  Discrete  Wavelet  Transform 


Discrete  wavelet  transform  was  implemented  as  a  perfect  reconstruction  filter  bank  [54].  In 
two  dimensions,  the  transform  is  obtained  by  applying  the  one-dimensional  transform 
separately  along  each  dimension.  Level  m  G  IV  of  the  transformed  image  matrix  s(nx,ny) 
is  therefore 


Uy)  ((>Am—  l5(n.3;,  Tly)  *  dfox))  4-2  *  9{Ply)')v 2) 
^y)  =  (O^m— l5(n.3;,  Tly)  *  g{nx))^2  *  h(jly)) _|_2, 

^y)  =  ((A-m— 15(72.3;,  Tly)  *  h(nx))±2  *  ^ (^1/) )j-2 s  &nd 
•Am 5 (n.3;,  Tly)  =  ((A771— 15(72.3;,  Tly)  ♦  h(nx))^2  *  h(jly))\2) 


where 


*A-qS{tIX)  Tly)  —  5(72.3;,  72y) 7 


(98) 


{yVis(nx,ny),W^s(nx,ny),W^s(rix,ny)}rne[1M]  and  AMs(nx,ny)  are  detail  and 
approximation  coefficients  for  M  levels  of  analysis,  respectively,  and  g(n)  and  h(n)  are  the 
decomposition  filters. 
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Reconstruction  of  the  approximation  coefficients  at  the  previous  level  is  given  by 

Am-is(nx,  ny )  =  ((W^s(nx,  ny))^2  *  g(ny )  +  (W*  s(nx,  ny))t2  *  h(ny))f2  *  g(nx)+ 
+((Wms(nx,  2  *  g(ny)  -f-  ( ~A.ms{nx ,  nyy)^2  *  h(nyy)^2  * 

with  g(n)  and  h(n)  being  the  reconstruction  filters. 

Two-dimensional  wavelets  associated  with  separable  filter  banks,  such  as  the  one  just 
described,  were  constructed  from  tensor  products  of  two  one-dimensional  multiresolution 
analyses  (wavelets  are  products  of  one-dimensional  wavelets  and  scaling  functions)  [17]. 
Note  also  that,  due  to  oversimplified  initialization  (98),  the  discrete  wavelet  transform  may 
be  a  pretty  poor  approximation  to  samples  of  the  continuous  wavelet  transform. 

We  will  limit  ourselves  to  FIR  filters  (i.e.,  compactly  supported  wavelets).  In  our 
experiments,  we  will  use  orthogonal  wavelet  transform  with  DAUB12  wavelet  [17],  and 
biorthogonal  wavelet  transform  with  Bior6.8  wavelet  from  MATLAB  Wavelet  Toolbox. 

2.3.3  Comparison  of  Transforms 

Image  fusion  is  performed  in  the  transform  space  by  computing  local  statistics  across  the 
decomposition  scales,  and  reconstructing  from  fused  transform  coefficients.  Typical  size  of 
neighborhood  is  between  a  single  pixel  and  5  by  5  area  with  some  loss  of  contrast  reported 
for  the  latter  [49].  In  general,  the  chosen  size  of  the  area  represents  a  tradeoff  between 
sharpness  and  introduction  of  artifacts.  In  order  to  achieve  our  goal  of  maximum  contrast 
with  minimum  artifacts,  we  limit  the  neighborhood  to  a  single  pixel  (maximizing  contrast) 
and  try  to  choose  the  most  appropriate  transform  (minimizing  artifacts). 

After  the  transform  decompositions  of  images  to  be  fused  is  performed,  the  corresponding 
transform  coefficients  of  the  images  are  combined  according  to  the  fusion  rule  into  a  new 
set  of  transform  coefficients  from  which  the  fused  result  is  reconstructed.  As  a  fusion  rule, 
we  used  the  maximum  magnitude  rule  (at  each  position  and  scale  of  the  transforms  the 
coefficient  with  greatest  magnitude  is  chosen)  for  the  gradient  pyramid,  the  orthogonal 
wavelet  transform  (DAUB12),  and  the  biorthogonal  wavelet  transform  (Bior6.8),  and  the 
maximum  oriented  energy  rule  (at  each  position  and  scale  of  the  transforms  the  coefficient 
with  greatest  local  oriented  energy  is  selected)  for  steerable  wavelet  transform  [46]. 

Our  first  two  experiments  used  the  phantom  shown  in  Figure  13.  Image  matrix  has 
dimensions  512  by  512,  and  fusion  was  performed  between  the  original  and  shifted 
phantom. 

First,  the  image  to  be  fused  with  the  one  from  Figure  13  was  generated  by  shifting  the 
phantom  one  sample  vertically  towards  the  top  of  the  image.  The  ideal  result  of  fusion 
should  contain  no  double  lines  or  other  artifacts  (the  distance  between  the  corresponding 
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Figure  13:  Phantom  used  for  comparisons  of  different  transforms  for  image  fusion. 
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Table  6:  Performance  of  fusion  algorithms  based  on  four  different  transforms. 


Transform 

MSE 

MAE 

Gradient  pyramid 

28.66 

13 

Orthogonal  WT 

0.14 

8 

Biorthogonal  WT 

0.21 

7 

Steerable  WT 

0.05 

6 

points  in  two  images  is  one  pixel,  so  that  the  algorithm  should  merge  shifted  features  into  a 
single  feature) .  Figure  14  shows  the  phantom  with  the  surrounding  area  in  the  fused 
images.  The  resulting  images  from  all  four  transforms  were  clipped  (pixel  values  above  the 
upper  limit  of  the  gray  level  range  were  mapped  to  white)  rather  than  scaled.  Please  note 
the  artifacts  present  when  the  orthogonal  and  biorthogonal  wavelet  transforms  were  used. 
The  latter  produced  a  slightly  better  result,  although  artifacts  due  to  aliasing  and  tensor 
product  representation  made  both  fused  images  unacceptable. 

Our  second  experiment  differed  from  the  first  one  only  in  the  fact  that  the  phantom  from 
Figure  13  was  shifted  by  five  samples.  Here,  the  shift  is  large  enough  that  features  from 
both  images  may  be  present  in  the  fused  image.  The  phantom  with  its  surroundings  in  the 
results  of  fusion  using  the  four  transforms  is  demonstrated  in  Figure  15.  Please  note  that 
the  situation  is  similar  to  the  one  in  the  first  experiment.  Both  orthogonal  and 
biorthogonal  wavelet  transforms  exhibited  obvious  artifacts,  while  the  redundant  gradient 
pyramid  and  steerable  dyadic  wavelet  transform  performed  well. 

The  third  experiment  was  geared  towards  a  quantitative  comparison  of  the  four  transforms 
for  image  fusion.  Similar  to  [52],  we  blur  different  parts  of  the  image  and  then  fuse  them  in 
such  a  way  that  each  blurred  part  is  fused  with  its  original  counterpart.  The  ideal  result  of 
fusion  would  be  the  original  image.  Figure  16  shows  the  original  512  by  512  mammogram, 
100  by  100  area  of  interest,  and  two  blurred  images  to  be  fused.  Mean-square  error  (MSE) 
and  maximum  absolute  error  (MAE)  between  the  result  of  fusion  and  the  original  image 
were  computed  for  all  four  methods.  Table  6  summarizes  the  results.  Let  us  remark  that 
these  results  are  rather  typical;  on  a  variety  of  images,  wavelet  based  methods  were  very 
close,  while  consistently  outperforming  the  gradient  pyramid  according  to  the  two  criteria. 
No  significant  artifacts  were  noticed  in  the  “blurring  experiments.”  By  comparing  the 
extracted  regions  with  the  original  ones,  we  subjectively  rated  the  transforms  used  for 
fusion  as:  (1)  steerable  dyadic  wavelet  transform,  (2)  biorthogonal  wavelet  transform,  (3) 
orthogonal  wavelet  transform,  and  (4)  gradient  pyramid.  Again,  the  differences  between 
different  types  of  wavelet  transforms  were  minor,  whereas  the  gradient  pyramid  lagged 
behind.  Regions  of  interest  corresponding  to  the  one  from  Figure  16(b)  are  shown  in  Figure 
17.  Wavelet  based  methods  produced  results  that  are  visually  close  to  the  original  region, 
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(CJ  W 


Figure  14:  Image  fusion  of  phantoms  shifted  by  one  sample,  (a)  Gradient  pyramid,  (b) 
Orthogonal  wavelet  transform,  (c)  Biorthogonal  wavelet  transform,  (d)  Steerable  dyadic 
wavelet  transform. 
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(c)  (d) 

Figure  15:  Image  fusion  of  phantoms  shifted  by  five  samples,  (a)  Gradient  pyramid,  (b) 
Orthogonal  wavelet  transform,  (c)  Biorthogonal  wavelet  transform,  (d)  Steerable  dyadic 
wavelet  transform. 
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(c)  (d) 


Figure  16:  (a)  Mammogram  with  area  of  interest  delineated,  (b)  Area  of  interest,  (c)  Image 
from  (a)  with  left  half  blurred,  (d)  Image  from  (a)  with  right  half  blurred. 
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gradient  pyramid,  however,  caused  loss  of  sharpness. 
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(a) 


(b) 


(c)  (d) 

Figure  17:  Region  of  interest  corresponding  to  the  one  from  Figure  16(b)  extracted  from  the 
fused  images  using:  (a)  gradient  pyramid,  (b)  orthogonal  wavelet  transform,  (c)  biorthogonal 
wavelet  transform,  and  (d)  steerable  dyadic  wavelet  transform. 
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3  Conclusions 


During  the  first  year,  we  have  made  significant  progress  in  the  development  of  a 
methodology  for  improving  the  mammographic  viewing  environment  by  steerable 
multiscale  transforms.  We  constructed  a  new  transform  that  does  not  introduce  artifacts 
due  to  translation  and  rotation  invariance,  which  are  inherent  to  traditional  wavelet 
analyses,  and  demonstrated  its  usefulness  for  image  fusion. 

We  extended  the  one-dimensional  discrete  dyadic  wavelet  transform  to  higher-order 
derivatives  and  even-order  spline  functions  and  developed  an  improved  initialization 
procedure.  Comparison  to  the  originally  employed  initialization  [11]  showed  significantly 
better  performance  of  our  procedure  for  finer  scales  of  analysis. 

We  developed  several  two-dimensional  transforms.  All  of  them  were  derived  with  the  goal 
of  eliminating  artifacts  due  to  lack  of  translation  and  rotation  invariance.  We  presented  a 
two-dimensional  discrete  dyadic  wavelet  transform  with  a  first-derivative  wavelet  as  an 
extension  of  the  one  originally  proposed  in  [11],  a  two-dimensional  discrete  dyadic  wavelet 
transform  with  a  second-derivative  wavelet  that  can  approximate  Laplacian  of  a  Gaussian, 
and  a  steerable  dyadic  wavelet  transform  implemented  in  a  near-perfect  reconstruction 
filter  bank  which  may  be  preferred  when  there  is  a  need  for  orientation  processing  along 
equally  spaced  angles.  We  derived  a  multiscale  spline  derivative-based  transform  which 
uses  x-y  separable  filters  in  a  perfect  reconstruction  filter  bank  and  enables  fast  translation 
and  rotation-invariant  directional  analysis  of  images. 

We  compared  the  steerable  dyadic  wavelet  transform  with  the  gradient  pyramid, 
orthogonal  wavelet  transform,  and  biorthogonal  wavelet  transform  for  fusion  of  features 
relevant  to  mammography.  During  our  experiments,  the  steerable  dyadic  wavelet  transform 
was  exhibiting  a  combination  of  best  properties  of  the  gradient  pyramid  and  of  the 
orthogonal/biorthogonal  wavelet  transform.  The  steerable  dyadic  wavelet  transform 
behaved  similarly  as  the  gradient  pyramid  in  a  sense  that  it  did  not  introduce  artifacts 
commonly  present  when  the  orthogonal  and  biorthogonal  wavelet  transforms  were  used.  At 
the  same  time,  the  steerable  dyadic  wavelet  transform  outperformed  the  gradient  pyramid 
by  acting  like  the  orthogonal  and  biorthogonal  wavelet  transforms  in  terms  of 
mathematical  criteria  and  sharpness  in  the  fused  image. 

Our  future  work  includes  building  of  a  unified  framework  for  creation  of  a  superior  local 
mammographic  viewing  environment  via  steerable  dyadic  wavelet  transform,  and  for  fusion 
of  locally  obtained  features  for  a  better  global  mammographic  viewing  environment. 
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